Регуляризация обратной задачи бигармонического уравнения
Список литературы
1. Тихонов А.Н., Самарский А.А. «Уравнения математической физики». -Главная редакция физико-математической литературы издательства «Наука», - М., 1966.
2. Тихонов А.Н., Арсенин В.Я. «Методы решения некорректных задач». - Главная редакция физико-математической литературы издательства «Наука», - М., 1974.
3. Тихонов А.Н., Гончарский А.В., Степанов В.В., Ягола А.Г. «Численные методы решения некорректных задач». - М., Наука, 1990.
4. Воеводин В.В., Кузнецов Ю.А. «Матрицы и вычисления». - М.: Наука, 1984г.
5. Вержбицкий В.М. «Численные методы. Математический анализ и обыкновенные дифференциальные уравнения». - М.: Высшая школа, 2001г.
6. Верлань А.Ф., Сизиков В.С. «Интегральные уравнения: методы, алгоритмы, программы», - Киев, 1986
Приложение
Описание переменных в программе:
n - число точек деления.
K1, K2, К3, К4 - функция ядер.
F - правая часть.
UTkm - точное решение,
UR - регуляризованное решение.
DU - невязка.
Текст программы:
clear;clc;disp("INVERSE PROBLEM for BIGARMONIC EQUATION:");
disp("d^4U/dx^4+2d^4U/(dxdy)^2+d^4U/dy^4=0 for(x,y)in CR");
disp("U(x,y)=FI(t),dU/dn=PCI(t) on L in CR CALCULATE U(x,y)");
disp("SYSTEM 1D-FR1:intg{s0,s1;K(t,s)z(s)ds}=FP(t),t0<=t<=t1");
disp("K(t,s)=[K1(t,s) K2(t,s); K3(t,s) K4(t,s)]");
disp("FP(t)=[FI(t);PSI(t)]");
N=input("ENTER N-NUMBER of KNOTS for FR1");
delta=input("Enter delta-error for FI(t),PCI(t) on L");
eps=input("Enter eps(<0.0001)for nev<=eps");
disp("R-radius for CR");
R=input(Enter R(test R=4));
disp ("(xp,yp)-center L as xp^2+yp^2<R^2");
xp=input (Enter xp);
yp=input (Enter yp);
RL2=(R-sqrt (xp^2+yp^2))^2
disp("Enter a1,b1-parameter for L as a1^2+b1^2<RL2");
a1=input(Enter a1-(test a1=3 for xp,yp=0));
b1=input(Enter b1-(test b1=2 for xp,yp=0));
s0=0; s1=2*%pi; t0=0; t1=2*%pi; N1=N-1;
//--------------------------------------------------------------
//disp("PRAM");hs=(s1-s0)/N;ht=(t1-t0)/N;for j=1:N;A(j)=hs;end;P=input(ENTER P={0;0.5;1});
disp("TRAP");P=1;hs=(s1-s0)/N1;ht=(t1-t0)/N1;A(1)=hs/2;A(N)=hs/2;for j=2:N1;A(j)=hs;end
//---------------Кривая L в CR ---------------------------------
function z=xL(t), z=a1*cos(t)+xp, endfunction
function z=yL(t), z=b1*sin(t)+yp, endfunction
//Производные параметров кривой L по x,y
function z=dxL(t), z=-a1*sin(t), endfunction
function z=dyL(t), z=b1*cos(t), endfunction
//Косинус(alfa) внешней нормали к кривой L
function z=cosal(t),z=dyL(t)/sqrt((dxL(t))^2+(dyL(t))^2),
endfunction
//Косинус(beta) внешней нормали к кривой L
function z=cosbe(t),z=-dxL(t)/sqrt((dxL(t))^2+(dyL(t))^2),
endfunction
//Функции для расчётных формул(из теории)
function z=a(x,y), z=x^2+y^2-R^2,endfunction
function z=d(x,y,s),
z=(x-R*cos(s))^2+(y-R*sin(s))^2,
endfunction
function z=c1(x,s), z=x-R*cos(s), endfunction
function z=c2(y,s), z=y-R*sin(s),endfunction
function z=b2(x,y,s), z=R-x*cos(s)-y*sin(s), endfunction
function z=K1(x,y,s),
z=-a(x,y)^2/(4*%pi*R*d(x,y,s)),
endfunction
function z=K2(x,y,s),
z=a(x,y)^2*b2(x,y,s)/(2*%pi*R*d(x,y,s)^2),
endfunction
//Производная K1(x,y,s) по х
function z=DK1X(x,y,s),
z=-(2*x*d(x,y,s)*a(x,y)-c1(x,s)*a(x,y)^2)/(2*%pi*d(x,y,s)^2*R),
endfunction
//Производная K1(x,y,s) по у
function z=DK1Y(x,y,s),
z=-(2*y*d(x,y,s)*a(x,y)-c2(y,s)*a(x,y)^2)/(2*%pi*d(x,y,s)^2*R),
endfunction
//Производная K1(x,y,s) по внешней нормали L
function z=K3(x,y,s,t),
z=DK1X(x,y,s)*cosal(t)+DK1Y(x,y,s)*cosbe(t), endfunction
//Производная K2(x,y,s) по х
function z=DK2X(x,y,s),
z=2*(d(x,y,s)*cos(s)+2*b2(x,y,s)*c1(x,s))*K1(x,y,s)/d(x,y,s)^2-
2*(b2(x,y,s)*DK1X(x,y,s))/d(x,y,s),
endfunction
//Производная K2(x,y,s) по у
function z=DK2Y(x,y,s),
z=2*(d(x,y,s)*sin(s)+2*b2(x,y,s)*c2(y,s))*K1(x,y,s)/d(x,y,s)^2-
2*(b2(x,y,s)*DK1Y(x,y,s))/d(x,y,s),
endfunction
//Производная K2(x,y,s) по внешней нормали L
function z=K4(x,y,s,t),
z=DK2X(x,y,s)*cosal(t)+DK2Y(x,y,s)*cosbe(t),
endfunction
//Точное решение Бигармонического уравнения
function z=UT(x,y), z=(x^2+y^2)/2, endfunction
//FI(t)-значение решения UT(x,y) на L
function z=FI(t), z=UT(xL(t),yL(t)),
endfunction
//Производная UT(x,y) по х
function z=UTx(x,y), z=x, endfunction
//Производная UT(x,y) по у
function z=UTy(x,y), z=y, endfunction
//PSI(t)-производная UT(x,y)по внешней нормали к L
function z=PSI(t),
z=UTx(xL(t),yL(t))*cosal(t)+UTy(xL(t),yL(t))*cosbe(t),
endfunction
//Построение С.Л.А.У.(дискретизация системы FR)
for i=1:N
t(i)=t0+(i-P)*ht;
for j=1:N
s(j)=s0+(j-P)*hs;
M1(i,j)=K1(xL(t(i)),yL(t(i)),s(j))*A(j);
M2(i,j)=K2(xL(t(i)),yL(t(i)),s(j))*A(j);
M3(i,j)=K3(xL(t(i)),yL(t(i)),s(j),t(i))*A(j);
M4(i,j)=K4(xL(t(i)),yL(t(i)),s(j),t(i))*A(j);
F1(i)=FI(t(i))+delta*(2*rand(t(i))-1);
F2(i)=PSI(t(i))+delta*(2*rand(t(i))-1);
end; end;
F=[F1;F2]; M=[M1 M2;M3 M4];
//Регуляризация Тихонова А.Н. С.Л.А.У.
K=10; q=0.1;
zp=(M*M+q^7*eye(2*N,2*N))M*F;
mu=norm(M*zp-F)^2;
for n=1:K;
alfa=q^n;
z=(M*M+alfa*eye(2*N,2*N))M*F;
nev=abs(norm(M*z-F)^2-delta^2-mu);
if nev<=eps then break
else
end; end;
disp("Exit alfa,nev when nev<=eps");
alfa=alfa
nev=nev
S=6; x0=-2; x1=2; y0=-2; y1=2;
for i=1:S
for j=1:S
xk(i)=x0+(i-1)*(x1-x0)/(S-1);
yk(j)=y0+(j-1)*(y1-y0)/(S-1);
Ut(i,j)=UT(xk(i),yk(j));
us1=0; us2=0;
for k=1:N
us1=us1+A(k)*z(k)/d(xk(i),yk(j),s(k));
us2=us2+A(k)*z(k+N)*b2(xk(i),yk(j),s(k))/d(xk(i),yk(j),s(k))^2;
end;
//Приближённое решение UR(x,y) уравнения
UR(i,j)=a(xk(i),yk(j))^2/(2*%pi*R)*((-0.5)*us1+us2);
end; end; S=S
S1=input(Enter 1<=S1<=S for rezults); m=1;
for i=1:S1
for j=1:S1
RES(m,1)=xk(i); RES(m,2)=yk(j);
RES(m,3)=Ut(i,j); RES(m,4)=UR(i,j);
RES(m,5)=UR(i,j)-Ut(i,j);
m=m+1; end; end; DUR=UR-Ut;
disp("EXIT REZULTS");
disp(" x y UT UR DUR ");
format(v,10); disp(RES); H=100;
//Графики(поверхности)
for j=1:H
tL(j)=t0+(j-1)*(t1-t0)/H;
XL(j)=xL(tL(j));
YL(j)=yL(tL(j));
CRx(j)=R*cos(tL(j));
CRy(j)=R*sin(tL(j));
end;
subplot(2,2,1)
plot3d(xk,yk,Ut)
xtitle("ТОЧНОЕ РЕШЕНИЕ")
subplot(2,2,2)
plot3d(xk,yk,UR)
xtitle("ПРИБЛИЖЁННОЕ РЕШЕНИЕ")
subplot(2,2,3)
plot3d(xk,yk,DUR)
xtitle("ПОГРЕШНОСТЬ")
subplot (2,2,4)
plot (XL,YL)
plot (CRx,CRy)
xtitle ("КРИВАЯ L в КРУГЕ CR")
Тестовый пример
В качестве замкнутой гладкой кривой L выберем
x = a*cos(t)
y = b*sin(t)
Рассмотрим однородную задачу:
Для иллюстрации результата выберем произвольные 9 точек.
u(х,у) =(x+y)/2.
Возьмем кривую с параметрами a=3, b=4; R=5; количество точек разбиения n=30: