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

求解波动方程数值解的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-a^2*u_xx=f(x,t) 0<x<1,0<t<T

%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);

本文来源：https://www.bwwdw.com/article/1yyj.html