НОУ ИНТУИТ | Моделирование систем. Лекция 9: Планирование и обработка результатов пассивного эксперимента

Учитесь и получайте официальные документы БЕСПЛАТНО. Вы можете поддержать наш проект.

Твой путь к знаниям!

Опубликован: 30.11.2010 | Уровень: специалист | Доступ: платный | ВУЗ: Мордовский государственный университет имени Н.П. Огарева

|

Вам нравится? Нравится 32 студентам

| Поделиться |

Поддержать программу

2. Аппроксимация функций степенными полиномами

С помощью метода наименьших квадратов возможно определять параметры уравнения кривых, заданных некоторыми своими точками, при условии, что вид уравнения известен, а неизвестными выступают параметры этого уравнения. При этом всякая аналитическая функция может быть приближена полиномом. Поэтому для подгонки нелинейных данных возможно использовать построение полинома методом наименьших квадратов. Для более точной подгонки степень полинома приходится увеличивать. Правда, если данные не проявляют полиномиальной природы, то в результате получается кривая, которая будет сильно осциллировать [15]. Этот феномен называется полиномиальное раскачивание [15, с. 301]. Оно явно наблюдается у полиномов высокой степени. Из этих соображений полиномы степени 6 или выше редко используются, если известно, что истинная функция, выражающая зависимость, является полиномом.

Но в то же время приближение данных полиномами служит основным инструментом для случая, когда уравнение кривой неизвестно совсем. Тогда в задачу входит определение коэффициентов (параметров) выбранного полинома, которым осуществляется подгонка. В любом случае метод наименьших квадратов дает достаточно приемлемые результаты. Более того, когда данные представляют собой некоторую "россыпь", с помощью метода наименьших квадратов возможно провести кривую между этими данными, которые наилучшим образом минимизируют квадрат невязки по этим данным.

Пример 2. Для функции вида $y=-5+1.2x+3.2x^3-12.4\sqrt{x}$ известны несколько значений точек на плоскости. По этим значениям рассчитайте параметры (коэффициенты) заданной функции.

Программный код решения примера:

clear,clc
try
   global hep
   delete(hep);
end
 
options.Resize = 'on';
options.WindowStyle = 'normal';
options.Interpreter = 'tex';
f = inputdlg({sprintf('\\bfВектор базисных функций:'),...
sprintf('\\bfОбласть определения:'),...
sprintf('\\bfЧисловые значения параметров b0,b1,...'),...
sprintf('\\bfЗначения аргумента (не менее 4 значений):...')},...
'ИСХОДНЫЕ ДАННЫЕ',1,...
{'x;x^3;sqrt(x)','[0 2]','-5;1.2;3.2;-12.4','0,0.12,0.7,1.34,0.45,1.678 '},options)

if isempty(f)
errordlg('Возможно, нажата клавиша ''Cansel'' ',...
            'Ошибка');
        return    
end
y1 = sym(char(f(1)));
x = str2num(char(f(2)));
b = str2num(char(f(3)));
es = str2num(char(f(4)));
  Y2 = char(y1);
q = find(Y2 == 'x');
if length(q) + 1 ~= length(str2num(char(f(3))) )
        errordlg('Число введенных параметров не соответствует числу параметров введенной функции',...
            'Ошибка ввода параметров');
        return
    end
%---------------- Общее число точек ---------------------
    if length(es) == 1
        if es < length(b)
             errordlg(sprintf('Общее число точек должно быть не менее числа %d - числа параметров',length(b)),...
                 'Ошибка ввода числа точек функции');
             return
         end
          fprintf('\n ОПРЕДЕЛЕНИЕ ПАРАМЕТРОВ ЗАДАННОЙ ФУНКЦИИ\n');
 fprintf('\t Вид параметрических функций:\n\t %s\n',Y2);
fprintf('\t Область определения функции: [%g %g]\n',x(1),x(2));
fprintf('\t Значения введенных параметров:\n');

for J = 1:length(b)
fprintf('\t%g\t',b(J)');
end
fprintf('\n');
disp('-----------------------------------------------------')
er1 = inline(sprintf('[1,%s]',char(f(1))));
erV1 = vectorize(inline(sprintf('[1,%s]',char(f(1)))));
erV = vectorize(inline(sprintf('[%s]',char(f(1)))));
con = char(erV);
n = find(con == ';');
con(n) = ',';
con1 = char(er1);
n1 = find(con1 == ';');
con1(n1) = ',';
symV = sym(con1);
wone = inline(con);
V =x(1) + (x(2) - x(1))*rand(es,1);
fprintf('\t Значения аргумента и параметрических функций:\n');
 
X = [ones(length(V),1),wone(V)];
disp([V,X])
XX = [];
for I = 1:length(V)
    XX = [XX;b'.*X(I,:)];
end
Y = sum(XX,2);
fprintf('\t Значения аргумента и заданной функции:\n');
disp([V,Y])
B = regress(Y,X);
fprintf('\t Расчетные значения параметров заданной функции:\n');
for J = 1:length(B)
    fprintf('\t b%d = %g\n',J - 1, B(J));
end
y = [];
for J = 1:length(symV)
    if J < length(symV)
        if J == 1
         y = [y,sprintf('b%d%s + ',J - 1,'')];
        else
    y = [y,sprintf('b%d*%s + ',J - 1,char(symV(J)))];
end
elseif J == length(symV)
    y = [y,sprintf('b%d*%s',J - 1,char(symV(end)))];
end
end
fprintf('\n\t Заданный вид функциональной зависимости:\n\t y = %s\n',y)
helpdlg(sprintf('Результаты выполнения программы смотрите в командном окне'),' ')
%----------------- Несколько точек ------------------------------
elseif length(es) >= length(b)
    if min(es) < x(1) | max(es) > x(2)
    errordlg(sprintf('Значения аргумента функции выходят за область определения [%g, %g]',x(1),x(2)),...
            'Ошибка ввода значений аргумента функции');
        return
    end
     fprintf('\n ОПРЕДЕЛЕНИЕ ПАРАМЕТРОВ ЗАДАННОЙ ФУНКЦИИ\n');
 fprintf(' Вид параметрических функций:\n\t %s\n',Y2);
fprintf(' Область определения функции: [%g %g]\n',x(1),x(2));
disp('-------------------------------------------------------------')
     er1 = inline(sprintf('[1,%s]',char(f(1))));
erV1 = vectorize(inline(sprintf('[1,%s]',char(f(1)))));
erV = vectorize(inline(sprintf('[%s]',char(f(1)))));
con = char(erV);
n = find(con == ';');
con(n) = ',';
con1 = char(er1);
n1 = find(con1 == ';');
con1(n1) = ',';
symV = sym(con1);
wone = inline(con);
V = es';
fprintf(' Значения аргумента и значения параметрических функций:\n');
X = [ones(length(V),1),wone(V)];
disp([V,X])
XX = [];
for I = 1:length(V)
    XX = [XX;b'.*X(I,:)];
end
Y = sum(XX,2);
fprintf(' Значения аргумента и значения заданной функции:\n');
disp([V,Y])
B = regress(Y,X);
fprintf(' Расчетные значения параметров заданной функции:\n');
for J = 1:length(B)
    fprintf('\t b%d = %g\n',J - 1,B(J));
end

y = [];
for J = 1:length(symV)
    if J < length(symV)
        if J == 1
         y = [y,sprintf('b%d%s + ',J - 1,'')];

        else
    y = [y,sprintf('b%d*%s + ',J - 1,char(symV(J)))];
end

elseif J == length(symV)
    y = [y,sprintf('b%d*%s',J - 1,char(symV(end)))];
end
end
fprintf('\n Заданный вид функциональной зависимости:\n\t y = %s\n',y);

global hep
hep = helpdlg(sprintf('Результаты выполнения программы смотрите в командном окне'), '');
else 
    errordlg(sprintf('Необходимо общее число аргументов не менее %d',length(b)), 'Ошибка!')
    return
end

Задание 2

Измените область определения параметрических функций, например, .
Постройте диаграмму заданной функциональной зависимости.
В программе задайте "свои" параметрические функции.
Определение параметров (коэффициентов) заданной функции произведите через решение нормального уравнения с помощью метода исключения Гаусса. Пример: .
В программе предусмотрите "зашумление" значений заданной функции по нормальному и равномерному вероятностным законам. Рассчитайте значения параметров (коэффициентов) заданной функции и подсчитайте относительную погрешность расчета коэффициентов. Подсчитайте также остаточную сумму квадратов;
Изучите функции , по справочной системе пакета MATLAB.

Пример 3. Функциональная экспериментальная зависимость задана массивом точек (значений). Произвести аппроксимацию данной функциональной зависимости степенным полиномом. Постройте графики функций.

Значения экспериментальной зависимости:

x	f(x)
0;	149.0;
5;	146.3;
10;	145.8;
15;	149.8;
20;	154.9;
25;	164.7;
30;	178.2;
35;	196.2;
40;	221.8;
45;	259.8;
50;	299.2;
55;	331.5;
60;	338.1;
65;	304.4;
70;	238.5;
72;	205.3;
74;	170.9;
75;	153.9;
76;	136.6;
78;	103.9;
80;	74.2;
82;	48.6;
84;	29.7;
85;	19.9;
86;	12.4;
88;	3.1;
90;	0;

Решение примера выполним с помощью стандартных функций MATLAB, таких как psi , finv , polyfit , polyval [8].

Программный код решения примера:

clear,clc,close all
x = ... [0;5;10;15;20;25;30;35;40;45;50;55;60;65;70;72;74;75;76;78;80;82;84;85;86;88;90];
f = ...
[149.0;146.3;145.8;149.8;154.9;164.7;178.2;196.2;221.8;259.8;299.2;331.5;338.1;304.4;238.5;205.3;170.9;153.9;
  136.6;103.9;74.2;48.6;29.7;19.9;12.4;3.1;0];   
n = length(x); %% length(x) == length(f)
 Lmin = f'*f; %% Квадратичный критерий - сумма квадратов
line(x,f,'linew', 2)
grid on
xlabel('\bf - - - - - - - - - x - - - - - - - - -');
ylabel('\bf f(x)');
 qprob = 0.4; % Уровень значимости
%%%% Выбор степени полинома
for m = 1 : n
    psi(:, m) = x.^(m-1); % psi - Специальная функция (полигамма)
    if m > 1
       for k = 1 : m-1
          psi(:, m) = psi(:, m) - ((psi(:,m))'*psi(:,k))*psi(:, k);
       end
    end
    psi(:, m) = psi(:,m)/norm(psi(:,m)); % нормировка
 
    bk = f'*psi(:,m);
    DOLD = Lmin/(n-m); %%%Старая дисперсия
    Lmin = Lmin - bk^2;
    Dnew = Lmin/(n-m-1); %%% Новая дисперсия
% fprintf('\n');
    if m > 1 % проверка степени начиная с 2
       if Dnew/DOLD < finv(qprob, n-m, n-m-1) 
% fprintf(' Степень %d следует учитывать\n', m-1);
           stp = m-1;
       else
% fprintf(' Степень %d не нужно учитывать\n', m-1);
 
if m > 2
 stp = m-2;
break;
else
  stp = m - 1;  
break;
end
        end
    end
end
 PO = polyfit(x, f, stp); 
fprintf('\n Коэффициенты аппроксимируещего полинома:\n'); 
for J = 1 : length(PO)
fprintf('%15g\n', PO(J));
end

fprintf('\n');
 yPO = polyval(PO, x);
 
line(x, yPO, 'marker', 'o', 'color','r')
title(sprintf('%s Аппроксимация экспериментальной кривой полиномом %d степени','\bf', stp));
legend('\bfЭксперимент','\bfАппроксимация', 'location', 'best')

Диаграммы заданной экспериментальной кривой и аппроксимирующего полинома показаны на рис. 9.2.

Рис. 9.2. Результат аппроксимации кривой полиномом

Задание 3

Проверьте результат выполнения программы при уменьшении и увеличении уровня значимости.
Постройте аппроксимирующую функцию непосредственно по аппроксимирующему полиному (без функции ).
Запишите в текстовый файл полученный аппроксимирующий полином. Имя файла задайте в виде compX.txt, где — номер компьютера, за которым выполняется лабораторная работа (1, 2, 3, ...).
Создайте массив значений функции на отрезке $[0; 4\pi]$ и произведите аппроксимацию полиномами следующих порядков: 2-го, 3-го, 4-го. Рассчитайте и постройте приведенную погрешность аппроксимации.

Примечание. В качестве приведенной погрешности примите разность в абсо-лютных значениях между аппроксимирующей функцией и истинной функ-ции, отнесенной к максимуму модуля истинной функции (например, для функции f(x) = cos(x) ).

Дальше >>

Моделирование систем

Моделирование систем

Планирование и обработка результатов пассивного эксперимента

2. Аппроксимация функций степенными полиномами

Вопросы и ответы

Студенты

Авторизоваться

Моделирование систем

Моделирование систем

Планирование и обработка результатов пассивного эксперимента

2. Аппроксимация функций степенными полиномами

Вопросы и ответы

Студенты