Наименьшее квадратное приближение на прямоугольнике [a, b] x [c, d]. Используйте функции 1, x, y, sin (x) и sin (y) в качестве основы - PullRequest
0 голосов
/ 13 января 2019

Мне нужно вычислить кофоцентричную матрицу b, но, независимо от того, какую функцию я пытаюсь использовать, созданная матрица A не имеет обратной, и поэтому A \ H не может быть получена.

Я проверил несколько раз, но я верю, что получил правильную матрицу правильно

function b = MeanSquareApprox(f,n,a,b,c,d)
%[a,b]x[c,d] - intervals
%n-number of nodes(each dimension)
%f-orginal function
h1=(b-a)/n;
h2=(d-c)/n;
x = a:h1:b; 
y = c:h2:d; 
exact=zeros(1,length(x));
for i=1:length(x)
    exact(1,i)=f(x(i),y(i));
end
[X,Y]=meshgrid(x,y);
F=f(X,Y);
surf(X,Y,F);

s1=sum(x);
s2=sum(y);
s3=sum(sin(x));
s4=sum(sin(y));

s5=sum(x.^2);
s6=sum(y.*x);
s7=sum(sin(x).*x);
s8=sum(sin(y).*x);

s9=sum(y.^2);
s10=sum(sin(x).*y);
s11=sum(sin(y).*y);

s12= sum(cos(x));
s13=sum(cos(x).^x);
s14=sum(cos(x).^y);
s15=sum(sin(x).*cos(x));
s16=sum(sin(y).*cos(x));

s17=sum(cos(y));
s18=sum(cos(y).*x);
s19=sum(cos(y).*y);
s20=sum(cos(y).*sin(x));
s21=sum(cos(y).*sin(y));

h1=sum(exact(1,:));
h2=sum(x.*exact(1,:));
h3=sum(y.*exact(1,:));
h4=sum(cos(x).*exact(1,:));
h5=sum(cos(y).*sin(y));

A=[length(x),s1,s2,s3,s4;s1,s5,s6,s7,s8;s2,s6,s9,s10,s11;s12,s13,s14,s15,s16;s17,s18,s19,s20,s21];
H=[h1;h2;h3;h4;h5];

disp(A);
disp(H);
b=A\H;
disp(b);
%g=@(x,y)(b(1)+b(2)*x+b(3)*y+b(4)*sin(x)+b(5)*sin(y))
end

enter image description here

1 Ответ

0 голосов
/ 14 января 2019

Вы можете использовать оператор Левенберга-Марквардта для улучшения стабильности деления матрицы.

Вместо

b = A\H;

Вы используете

b = (A+.001.*eye(size(A,1)))\H;

Параметр mu может быть установлен как очень маленькое значение, и все, что он делает, это стабилизирует деление матрицы. Для системы, которую вы показали ранее со значением mu 0,001, я получил:

    -0.104571468432964
    -0.414935581633317
    -0.414935563888507
      41460.5762052549
     -41457.2679510469

Есть несколько других способов сделать это, но это очень сложная тема, и это самый простой.

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...