(完整版)求解波动方程数值解的matlab程序隐式格式2021.docx

PAGE 1Word版本 PAGE 1 (完整版)求解波动方程数值解的matlab程序隐式格式2021 求解波动方程数值解的matlab程序隐式格式2021-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-a^2*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-1  U(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)));  end  for j=1:n  U(j,1)=feval(ft1,t(j));  U(j,m)=feval(ft2,t(j));  end  A=-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;  for i=3:n  for j=2:m-1  UU(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));  end  f1(1)=f1(1)+0.5*s^2*U(i,1);  f1(end)=f1(end)+0.5*s^2*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);