erweipaowu.m

在matlab环境下的二维双曲型和抛物型偏微分方程的数值解法

function [s1,s2,s3]=erweipaowu(f,a,b,h,t,s0)
%s1,s2,s3分别是数值解，精确解和误差在指定节点上的值，f,a,b是与方程有关的函数，h,t分别是空间和时间步长，s0则是精确解
r=t/(h*h);
m=1/h;
n=1/t;
for i=1:m+1
    x(i)=(i-1)*h;
    y(i)=(i-1)*h;
end
for k=1:n+1
    z(k)=(k-1)*t;
end
for i=2:m
    for j=2:m
        u(i,j,1)=a(x(i),y(j));
    end
end
for k=1:n+1
    for j=1:m+1
        u(1,j,k)=b(x(1),y(j),z(k));
        u(m+1,j,k)=b(x(m+1),y(j),z(k));
    end
    for i=2:m
        u(i,1,k)=b(x(i),y(1),z(k));
        u(i,m+1,k)=b(x(i),y(m+1),z(k));
    end
end
w1=ones(1,m-1);w2=ones(1,m-2);
A1=diag(w1);A2=diag(w2,-1);A3=diag(w2,1);
A=(1+r)*A1-r/2*A2-r/2*A3;
B=r*(1-r)/2*A1+r*r/4*A2+r*r/4*A3;
C=(1-r)*(1-r)*A1+r*(1-r)/2*A2+r*(1-r)/2*A3;
for k=1:n
    for j=2:m
        v1(1,j-1)=u(1,j,k+1)-r/2*(u(1,j-1,k+1)-2*u(1,j,k+1)+u(1,j+1,k+1));
        v1(m+1,j-1)=u(m+1,j,k+1)-r/2*(u(m+1,j-1,k+1)-2*u(m+1,j,k+1)+u(m+1,j+1,k+1));
        for i=1:m-1
            d1(i)=u(i+1,j-1,k);
            d2(i)=u(i+1,j,k);
            d3(i)=u(i+1,j+1,k);
        end
        d4(1)=r/2*v1(1,j-1)+r*r/4*u(1,j-1,k)+r*(1-r)/2*u(1,j,k)+r*r/4*u(1,j+1,k)+t/2*(f(x(2),y(j),z(k))+f(x(2),y(j),z(k+1)));
        d4(m-1)=r/2*v1(m+1,j-1)+r*r/4*u(m+1,j-1,k)+r*(1-r)/2*u(m+1,j,k)+r*r/4*u(m+1,j+1,k)+t/2*(f(x(m),y(j),z(k))+f(x(m),y(j),z(k+1)));
        for i=2:m-2
            d4(i)=t/2*(f(x(i+1),y(j),z(k))+f(x(i+1),y(j),z(k+1)));
        end
        d=B*d1'+C*d2'+B*d3'+d4';
        v=chase(A,d,m-2);
        for i=2:m
            v1(i,j-1)=v(i-1);
        end
    end
    for i=2:m
        u(i,1,k+1)=b(x(i),y(1),z(k+1));
        u(i,m+1,k+1)=b(x(i),y(m+1),z(k+1));
        d(1)=v1(i,1)+r/2*u(i,1,k+1);
        d(m-1)=v1(i,m-1)+r/2*u(i,m+1,k+1);
        for j=2:m-2
            d(j)=v1(i,j);
        end
        v=chase(A,d,m-2);
        for j=2:m
            u(i,j,k+1)=v(j-1);
        end
    end
end
for i=1:m+1
    for j=1:m+1
        for k=1:n+1
            u0(i,j,k)=s0(x(i),y(j),z(k));
        end
    end
end
E=u-u0;
s1=single(u(11,11,3:2:21));
s2=single(u0(11,11,3:2:21));    
s3=single(abs(E(11,11,3:2:21)));