求解波動(dòng)方程數(shù)值解的matlab程序隱式格

1、求解波動(dòng)方程數(shù)值解的 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,

2、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)%求解下問(wèn)題%u_tt-aA2*u_xx=f(x,t) 0<x<1,0<t<T%u(

3、0,t)=ft1,u(1,t)=ft2,%u(x,0)=fx1(x)%u_t(x,0)=fx2(x)%h離散x方向的步長(zhǎng)%k離散t方向的步長(zhǎng)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)+kA2/2*(aA2/hA2* .(feval(fx1,x(i+1)-2*feval(fx1,x(i)+feval(fx1,x(i-1)+feval

4、(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*sA2*ones(1,m-2);C=A;B=(1+sA2)*ones(1,m-2);UU=zeros(1,m-2);f1=UU;for i=3:nfor j=2:m-1UU(j-1)=f0(x(j),t(i);f1(j-1)=0.5*sA2*U(i-2,j-1)-(1+sA2)*U(i-2,j).+0.5*sA2*U(i-2,j+1)+2*U(i-1,j).+kA2*feval

5、(f,x(j),t(i-1);endf1(1)=f1(1)+0.5*sA2*U(i,1);f1(end)=f1(end)+0.5*sA2*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-1B1(i)=C(i)/(B(i)-A(i)*B1(i-1);

6、endY(1)=f(1)/B(1);for i=2:nY(i)=(f(i)-A(i)*Y(i-1)/(B(i)-A(i)*B1(i-1);endx1(n)=Y(n);for i=n-1:-1:1x1(i)=Y(i)-B1(i)*x1(i+1);endx=x1;function z=f0(x,t)%精確解函數(shù)z=exp(x+t);O(D拉格朗日插值的Matlab實(shí)現(xiàn)0 Matlab中沒(méi)有現(xiàn)成的拉格朗日插值函數(shù)，必須編寫(xiě)一個(gè)文件實(shí)現(xiàn)拉格朗日插值。設(shè)n個(gè)節(jié)點(diǎn)數(shù)據(jù)以數(shù)組 x0, y0輸入(注意 Matlab的數(shù)組下標(biāo)從1開(kāi)始)，m個(gè)插值點(diǎn)以數(shù) 組x輸入，輸出數(shù)組y為m個(gè)插值。編寫(xiě)一個(gè)名為lagrange.m的M文件：function y=lagrange(x0,y0,x);n=length(x0);m=length(x);for i=1:mz=