Question

求解二维抛物线型偏微分方程matlab程序

方程如下：

求h,
需要编出求解的matlab程序，谢谢！

Answer 1

function [u,x,y,t] = TDE(A,D,T,ixy0,bxyt,Mx,My,N)
%解方程 u_t = c(u_xx + u_yy) for D(1) <= x <= D(2), D(3)<= y <= D(4), 0 <= t<= T
%初值: u(x,y,0) = ixy0(x,y)
% 边界条件: u(x,y,t) = bxyt(x,y,t) for (x,y)cB
% Mx/My ： x轴和y轴的等分段数
% N ： t 轴的等分段数
dx = (D(2) - D(1))/Mx; x = D(1)+[0:Mx]*dx;
dy = (D(4) - D(3))/My; y = D(3)+[0:My]'*dy;
dt = T/N; t = [0:N]*dt;
%初始化u
for i = 1:Mx + 1
for j = 1:My +1
u(i,j) =ixy0(x(i),y(j));
end
end
rx = A*dt/(dx*dx);
rx1 = 1 + 2*rx; rx2 = 1 - 2*rx;
ry = A*dt/(dy*dy);
ry1 = 1 + 2*ry; ry2 = 1 - 2*ry;
for i = 1:Mx - 1 %(11.2.21a)
P(i,i) = ry1;
if i > 1
P(i - 1,i)= -ry; P(i,i-1) = -ry;
end
end
for j = 1:My - 1 %(11.2.21b)
Q(j,j) = rx1;
if j > 1
Q(j- 1,j) =-rx; Q(j,j - 1) = -rx;
end
end

for k = 1:N
u_1 = u; t =k*dt;
for i =1:Mx + 1 %边界条件
u(i,1) = feval(bxyt,x(i),y(1),t);
u(i,My+1) = feval(bxyt,x(i),y(My+1),t);
end
for j =1:My + 1
u(1,j)= feval(bxyt,x(1),y(j),t);
u(Mx+1,j) = feval(bxyt,x(Mx+1),y(j),t);
end
if mod(k,2)== 0
for i =2:Mx
j =2:My;
bx=[ry*u(i,1) zeros(1,Mx - 3) ry*u(i,My + 1)] +rx*(u_1(i-1,j)+ u_1(i + 1,j)) +rx2*u_1(i,j);
u(i,j) = linsolve(P,bx'); %(11.2.21a)
end
else
for j = 2:My
i =2:Mx;
by=[rx*u(1,j); zeros(My-3,1); rx*u(Mx + 1,j)] + ry*(u_1(i,j-1) + u_1(i,j + 1)) +ry2*u_1(i,j);
u(i,j) = linsolve(Q,by); %(11.2.21b)
end
end
end
u=u';