求解波动方程数值解的matlab程序隐式格式2010.doc

求解波动方程数值解的matlab程序 隐式格式2010-04-19 13:45function varargout=liu(varargin)a=1;T=1;a=1;b=0.5;h=1/20;k=1/40;f=inline(0,x,t);fx1=inline(exp(x);fx2=inline(exp(x);ft1=inline(exp(t);ft2=inline(exp(1+t);X,Y,Z=chfenmethed(f,fx1,fx2,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,fx1,fx2,ft1,ft2,a,T,h,k)%求解下问题%u_tt-a2*u_xx=f(x,t) 0x1,0tT%u(0,t)=ft1,u(1,t)=ft2,%u(x,0)=fx1(x)%u_t(x,0)=fx2(x)%h离散x方向的步长%k离散t方向的步长x=0:h:a;t=0:k:T;m=length(x);n=length(t);s=a*k/h;X,T=meshgrid(x,t);Z=zeros(n,m);U=zeros(n,m);for i=2:m-1U(1,i)=feval(fx1,x(i);U(2,i)=U(1,i)+k*feval(fx2,x(i)+k2/2*(a2/h2* .(feval(fx1,x(i+1)-2*feval(fx1,x(i)+feval(fx1,x(i-1)+feval(f,x(i),0);Z(2,i)=abs(U(2,i)-f0(x(i),t(2);endfor j=1:nU(j,1)=feval(ft1,t(j);U(j,m)=feval(ft2,t(j);endA=-0.5*s2*ones(1,m-2);C=A;B=(1+s2)*ones(1,m-2);UU=zeros(1,m-2);f1=UU;for i=3:n for j=2:m-1 UU(j-1)=f0(x(j),t(i); f1(j-1)=0.5*s2*U(i-2,j-1)-(1+s2)*U(i-2,j). +0.5*s2*U(i-2,j+1)+2*U(i-1,j). +k2*feval(f,x(j),t(i-1); end f1(1)=f1(1)+0.5*s2*U(i,1); f1(end)=f1(end)+0.5*s2*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);endfunction x=zgf(A,B,C,f)%解 b1 c1% a2 b2 c2% . . . *x=f% % an bnn=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);endY(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);endx1(n)=Y(n);for i=n-1:-1:1 x1(i)=Y(i)-B1(i)*x1(i+1); endx=x1;function z=f0(x,t)%精确解函数z=exp(x+t);