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

《求解波动方程数值解的matlab程序隐式格式2010.doc》由会员上传分享，免费在线阅读，更多相关内容在工程资料-天天文库。

1、求解波动方程数值解的matlab程序隐式格式2010-04-1913:45functionvarargout=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);shadingflat;xlabel('X',

2、'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-a^2*u_xx=f(x,t)0

3、=a*k/h;[X,T]=meshgrid(x,t);Z=zeros(n,m);U=zeros(n,m);fori=2:m-1U(1,i)=feval(fx1,x(i));U(2,i)=U(1,i)+k*feval(fx2,x(i))+k^2/2*(a^2/h^2*...(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)));endforj=1:nU(j,1)=feval(ft1,t(j));U(j,m)=feval(f

4、t2,t(j));endA=-0.5*s^2*ones(1,m-2);C=A;B=(1+s^2)*ones(1,m-2);UU=zeros(1,m-2);f1=UU;fori=3:nforj=2:m-1UU(j-1)=f0(x(j),t(i));f1(j-1)=0.5*s^2*U(i-2,j-1)-(1+s^2)*U(i-2,j)...+0.5*s^2*U(i-2,j+1)+2*U(i-1,j)...+k^2*feval(f,x(j),t(i-1));endf1(1)=f1(1)+0.5*s^2*U(i,1);f1(end)=f1(end)+0.5*s^2*U(i,m);U

5、(i,2:m-1)=zgf(A,B,C,f1);Z(i,2:m-1)=abs(U(i,2:m-1)-UU);endfunctionx=zgf(A,B,C,f)%解[b1c1%a2b2c2%...*x=f%%anbn]n=length(B);B1=zeros(1,n-1);Y=zeros(1,n);x1=zeros(1,n);B1(1)=C(1)/B(1);fori=2:n-1B1(i)=C(i)/(B(i)-A(i)*B1(i-1));endY(1)=f(1)/B(1);fori=2:nY(i)=(f(i)-A(i)*Y(i-1))/(B(i)-A(i)*B1(i-1));

6、endx1(n)=Y(n);fori=n-1:-1:1x1(i)=Y(i)-B1(i)*x1(i+1);endx=x1;functionz=f0(x,t)%精确解函数z=exp(x+t);