xqeuler.m

在matlab环境下的向前和向后euler算法价偏微分方程问题

function [E3]=xqeuler(t,h,f,z,a,b,c,s)
%h是空间步长，t是时间步长，f是方程右端项，z是u的初值函数，a、b分别是u的边界函数，c是对空间二阶偏导前系数的相反数,s是精确解函数。
m=1/h;n=1/t;
r=(c*t)/(h*h);
for i=1:m+1
    x(i)=(i-1)*h;
    u(i,1)=z(x(i));
end
for i=1:n+1
    y(i)=(i-1)*t;
    u(1,i)=a(y(i));
    u(m+1,i)=b(y(i));
end
A(1,1)=1-2*r;A(1,2)=r;A(m-1,m-2)=r;A(m-1,m-1)=1-2*r;
for i=2:m-2
    A(i,i)=1-2*r;A(i,i-1)=r;A(i,i+1)=r;
end
for k=1:n
    for i=1:m-1
        d(i)=u(i+1,k);
    end
    for i=2:m-2
        e(i)=t*f(x(i+1),y(k));
    end
    e(1)=t*f(x(2),y(k))+r*u(1,k);
    e(m-1)=t*f(x(m),y(k))+r*u(m+1,k);
    v=A*d'+e';
for i=2:m
    u(i,k+1)=v(i-1);
end
end
for i=1:m+1
    w(i)=s(x(i),1);
end
E3=abs(u(:,n+1)-w(:));