xheuler.m

来自「在matlab环境下的向前和向后euler算法价偏微分方程问题」· M 代码 · 共 31 行

function [E3]=xheuler(t,h,f,z,a,b,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=2:n+1
    for i=2:m-2
        d(i)=u(i+1,k-1)+t*f(x(i+1),y(k));
    end
    d(1)=u(2,k-1)+t*f(x(2),y(k))+r*u(1,k);
    d(m-1)=u(m,k-1)+t*f(x(m),y(k))+r*u(m+1,k);
    v=chase(A,d,m-2);
    for i=2:m
        u(i,k)=v(i-1);
    end
end
for i=1:m+1
    o(i)=s(x(i),1);
end
E3=abs(u(:,n+1)-o(:));