向后 Euler法求解抛物线方程的matlab程序

function varargout=liu(varargin) T=1;a=1;h=1/10;k=1/400; f=inline('0','x','t'); fx=inline('exp(x)'); ft1=inline('exp(t)'); ft2=inline('exp(1+t)');

[X,Y,Z]=chfenmethed(f,fx,ft1,ft2,a,T,h,k); mesh(X,Y,Z); shading flat;

xlabel('X','FontSize',14); ylabel('t','FontSize',14);

zlabel('error','FontSize',14); title('误差图');

function [X,T,Z]=chfenmethed(f,fx,ft1,ft2,a,T,h,k) %求解下问题

%u_t-a*u_xx=f(x,t) 0m=length(x); n=length(t); r=a*k/h^2;

[X,T]=meshgrid(x,t); Z=zeros(n,m); U=zeros(n,m); for i=1:m

U(1,i)=feval(fx,x(i)); end

for j=2:n

U(j,1)=feval(ft1,t(j)); U(j,m)=feval(ft2,t(j)); end

A=-r*ones(1,m-2);

B=(1+2*r)*ones(1,m-2); C=A;

UU=zeros(1,m-2); f1=UU;

for i=2:n for j=2:m-1

UU(j-1)=f0(x(j),t(i));

f1(j-1)=U(i-1,j)+k*feval(f,x(j),t(i)); end

f1(1)=f1(1)+r*U(i,1); f1(m-2)=f1(m-2)+r*U(i,m); U(i,2:m-1)=zgf(A,B,C,f1);

Z(i,2:m-1)=abs(U(i,2:m-1)-UU); end

function x=zgf(A,B,C,f) %解 [b1 c1

% a2 b2 c2

% . . . *x=f % % an bn] n=length(B);

B1=zeros(1,n-1); Y=zeros(1,n); x1=zeros(1,n); B1(1)=C(1)/B(1); for i=2:n-1

B1(i)=C(i)/(B(i)-A(i)*B1(i-1)); end

Y(1)=f(1)/B(1); for i=2:n

Y(i)=(f(i)-A(i)*Y(i-1))/(B(i)-A(i)*B1(i-1)); end

x1(n)=Y(n); for i=n-1:-1:1

x1(i)=Y(i)-B1(i)*x1(i+1); end x=x1;

function z=f0(x,t) %精确解函数 z=exp(x+t);