Применение кубатуры для регуляризации задачи идентификации значения оператора Лапласа

N=input('Enter N-number of knots for OX,OY');

delta=input('Enter delta-error for U(x,y)');

eps=input('Enter eps(<0.0001)for nev<=eps');

L=0;

 

//=========Концы отрезка для х,у и параметры А,B>0================

a=0;   b=1;     c=0;   d=2;    A=1; B=1; 

disp("PRAM");ht=(b-a)/N;hs=(d-c)/N;for j=1:N;At(j)=ht;As(j)=hs;end;P=input('Enter P={0;0.5;1}');

RU=2*rand(N,N)-ones(N,N);

 

//=========Ядро(функция Грина)для интегрирования по х=============

function w=G1(x,t),

if t<=x then    

w=(x-b)*(t-a)/(b-a);

elseif x<=t then

w=(x-a)*(t-b)/(b-a);

end;   endfunction

 

//=========Ядро(функция Грина)для интегрирования по у=============

function w=G2(y,s),

if s<=y then     

w=(y-d)*(s-c)/(d-c);

elseif y<=s then

w=(y-c)*(s-d)/(d-c);

end;   endfunction

 

//=========U(x,y),zt(x,y),d^2U(x,y)/dx^2,d^2U(x,y)/dy^2 в PDE=====

disp("Example"); deff('w=U(x,y)','w=x^3*y+y^4'); deff('w=d2Ux2(x,y)','w=6*x*y');deff('w=d2Uy2(x,y)','w=12*y^2'); deff('w=zt(x,y)','w=A*d2Ux2(x,y)+B*d2Uy2(x,y)+L*U(x,y)');

 

//=========Построение функции T(x,y)==============================

deff('w=T1(x,y)','w=((c-y)*d2Ux2(x,d)+(y-d)*d2Ux2(x,c))/(d-c)');  

deff('w=T2(x,y)','w=((a-x)*d2Uy2(b,y)+(x-b)*d2Uy2(a,y))/(b-a)');

deff ('w=T(x,y)','w=-A*T1(x,y)-B*T2(x,y)+L*U(x,y)');

 

//=========Построение функции f(x,y)==============================

deff('w=f1(x,y)','w=((a-x)*U(b,y)+(x-b)*U(a,y))/(b-a)');

deff('w=f2(x,y)','w=(c-y)*(U(x,d)+((a-x)*U(b,d)+(x-b)*U(a,d))/(b-a))/(d-c)');

deff('w=f3(x,y)','w=(y-d)*(U(x,c)+((a-x)*U(b,c)+(x-b)*U(a,c))/(b-a))/(d-c)');

deff('w=f(x,y)','w=f1(x,y)+f2(x,y)+f3(x,y)');

deff('w=W(x,y)','w=U(x,y)+f(x,y)');

 

//=========Численное интегрирование  для дискретизации W(x,y)======

hx=ht;  hy=hs;

for i=1:N;

for j=1:N;

x(i)=a+(i-P)*hx;       y(j)=c+(j-P)*hy;

end;   end;    

t=x;   s=y;   

N2=100;  

ht2=(b-a)/N2; hs2=(d-c)/N2;for j=1:N2;At2(j)=ht2;As2(j)=hs2;end;

for i=1:N;

for j=1:N;

S1=0;     S2=0;

for m=1:N2

t2(m)=a+(m-0.5)*ht2;      s2(m)=c+(m-0.5)*hs2; 

S1=S1+B*G1(x(i),t2(m))*W(t2(m),y(j))*At2(m);

S2=S2+A*G2(y(j),s2(m))*W(x(i),s2(m))*As2(m);

end;

u(i,j)=S1+S2+delta*RU(i,j); 

end;    end;

 

//=========Правая часть развёрнутой  системы=======================

for p=1:N;

for k=1:N;

q=(p-1)*N+k;

F(q)=u(p,k);

end;   end;

 

//=========4-индексная матрица G(p,k,i,j) в 2-индексную===========

for p=1:N;

for k=1:N;

q=(p-1)*N+k;

for i=1:N;

for j=1:N;

n=(j-1)*N+i;

 

//=========2-индексная матрица развёрнутой  системы по кубатуре====

Q(q,n)=G1(x(p),t(j))*G2(y(k),s(i))*At(j)*As(i);

end;   end;  end;    end;

 

//====Регуляризация С.Л.А.У. mu=inf{norm(Q*z-F)^2} all z in E(N^2)=

L=10; g=0.1;

zp=(g^8*eye(N^2,N^2)+Q'*Q)\(Q'*F);

mu=(norm(Q*zp-F))^2;

for m=1:L;

alfa=g^m;

z=(alfa*eye(N^2,N^2)+Q'*Q)\(Q'*F);

nev=abs(norm(Q*z-F)^2-delta^2-mu);

if nev<=eps then break

else

end;      end;

disp("Exit nev,alfa for nev<=eps");

nev=nev 

alfa=alfa

 

//=========Построение решений, погрешности  как матриц==============

for i=1:N;

  for k=1:N;

ZT(i,k)=zt(x(i),y(k)); 

H(i,k)=z(k+(i-1)*N);

ZR(i,k)=H(i,k)+T(x(i),y(k));

end;    end;  

DZR=ZR-ZT;

M=input('Enter 1<=M<=N for exit results'); m=1;     

for i=1:M

for j=1:M 

RES(m,1)=x(i);    RES(m,2)=y(j);    RES(m,3)=ZT(i,j);

RES(m,4)=ZR(i,j); RES(m,5)=DZR(i,j);m=m+1; end; end;

disp("   x        y       ZT        ZR      DZR    ");

format('v',7);            disp(RES);

 

//Графики

subplot(3,1,1)

plot3d(x,y,ZT)

xtitle("ZT-ТОЧНОЕ РЕШЕНИЕ,(x,y)-KNOTS")

subplot(3,1,2)

plot3d(x,y,ZR)

xtitle("ZR–РЕГУЛЯРИЗОВАННОЕ РЕШЕНИЕ,(x,y)-KNOTS")

subplot(3,1,3)

plot3d(x,y,DZR)

xtitle("ПОГРЕШНОСТЬ РЕГУЛЯРИЗАЦИИ,(x,y)-KNOTS")

 

 

Пример . Для функции при . Точное решение (для контроля) = . Погрешность задана в программе.

График 1

Здесь взято           Числовые результаты.

,                x        y       ZT           ZR        DZR

               0.      0.    0.             0.            0.

              0.      0.1    0.121      0.121       0.

и кубатуры  прямоугольников (р=1).           0.      0.2    0.496      0.496      0.

             0.      0.2    0.496     0.496    0.

              0.      0.3    1.161     1.161    0.

               0.      0.4    2.176     2.176      0.

              0.      0.5    3.625     3.625    0.

              0.      0.6    5.616     5.616    0.

               0.       0.7   8.281  8.281   0.

                 0.      0.8    11.776    11.776   0.

                 0.      0.9    16.281    16.281   0.

                 0.         1.         22.     22.        0.

                 0.      1.1    29.161    29.161  0.

                0.      1.2    38.016    38.016  0.