erweishuangqu.m

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

function [s1,s2,s3]=erweishuangqu(f,a,b,c,h,t,s0)
%s1,s2,s3分别是数值解，精确解和误差在指定节点上的值，f,a,b是与方程有关的函数，h,t分别是空间和时间步长，s0则是精确解
s=t/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=1:m+1
    for j=1:m+1
        u(i,j,1)=a(x(i),y(j));
        u(i,j,2)=a(x(i),y(j))+t*b(x(i),y(j))+(t*t)/2*f(x(i),y(j),0);
    end
end
for k=3:2:n+1
    for j=1:m+1
        u(1,j,k)=c(x(1),y(j),z(k));
        u(m+1,j,k)=c(x(m+1),y(j),z(k));
    end
    for i=2:m
        u(i,1,k)=c(x(i),y(1),z(k));
        u(i,m+1,k)=c(x(i),y(m+1),z(k));
    end
end
for k=4:2:n+1
    for j=2:m
        u(1,j,k)=c(x(1),y(j),z(k))+u(1,j,2)-c(x(1),y(j),z(2));
        u(m+1,j,k)=c(x(m+1),y(j),z(k))+u(m+1,j,2)-c(x(m+1),y(j),z(2));
    end
    for i=2:m
        u(i,1,k)=c(x(i),y(1),z(k))+u(i,1,2)-c(x(i),y(1),z(2));
        u(i,m+1,k)=c(x(i),y(m+1),z(k))+u(i,m+1,2)-c(x(i),y(m+1),z(2));
    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+s*s)*A1-s*s/2*A2-s*s/2*A3;
B=-4/(h*h)*A1+1/(h*h)*A2+1/(h*h)*A3;
for k=2:n
    for j=2:m
        v1(1,j-1)=-s*s/(2*t*t)*(u(1,j-1,k-1)+u(1,j+1,k-1)-2*u(1,j-1,k)-2*u(1,j+1,k)+u(1,j-1,k+1)+u(1,j+1,k+1))+(1+s*s)/(t*t)*(u(1,j,k-1)-2*u(1,j,k)+u(1,j,k+1));
        v1(m+1,j-1)=-s*s/(2*t*t)*(u(m+1,j-1,k-1)+u(m+1,j+1,k-1)-2*u(m+1,j-1,k)-2*u(m+1,j+1,k)+u(m+1,j-1,k+1)+u(m+1,j+1,k+1))+(1+s*s)/(t*t)*(u(m+1,j,k-1)-2*u(m+1,j,k)+u(m+1,j,k+1));
        for i=1:m-1
            d1(i)=u(i+1,j,k);
        end
        d2(1)=1/(h*h)*(u(2,j-1,k)+u(2,j+1,k)+u(1,j,k))+f(x(2),y(j),z(k))+s*s/2*v1(1,j-1);
        d2(m-1)=1/(h*h)*(u(m,j-1,k)+u(m,j+1,k)+u(m+1,j,k))+f(x(m),y(j),z(k))+s*s/2*v1(m,j-1);
        for i=2:m-2
            d2(i)=(u(i+1,j-1,k)+u(i+1,j+1,k))/(h*h)+f(x(i+1),y(j),z(k));
        end
        d=B*d1'+d2';
        v=chase(A,d,m-2);
        for i=2:m
            v1(i,j-1)=v(i-1);
        end
    end
    for i=2:m
        u1(i-1,1)=(u(i,1,k-1)-2*u(i,1,k)+u(i,1,k+1))/(t*t);
        u1(i-1,m+1)=(u(i,m+1,k-1)-2*u(i,m+1,k)+u(i,m+1,k+1))/(t*t);
        d(1)=v1(i,1)+(s*s)/2*u1(i-1,1);
        d(m-1)=v1(i,m-1)+(s*s)/2*u1(i-1,m+1);
        for j=2:m-2
            d(j)=v1(i,j);
        end
        v=chase(A,d,m-2);
        for j=2:m
            u1(i-1,j)=v(j-1);
        end
    end
    for i=2:m
        for j=2:m
            u(i,j,k+1)=(t*t)*u1(i-1,j)+2*u(i,j,k)-u(i,j,k-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(21,21,5:4:41));
s2=single(u0(21,21,5:4:41));    
s3=single(abs(E(21,21,5:4:41)));