tde.m

偏微分方程matlab m文件

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';