偏微分方程数值及matlab实验报告(8)

偏微分方程数值实验报告八实验题目：利用有限差分法求解.0)1(,0)1(),()()(uuxfxuxu真解为)1()(22xexux实现算法：对于两点边值问题,)(,)(,,d22buaulxfdxu(1)其中),(balfba),(为],[bal上的连续函数，，为给定常数.其相应的有限差分法的算法如下：1.对求解区域做网格剖分，得到计算网格.在这里我们对区间l均匀剖分n段，每个剖分单元的剖分步长记为nabh.2.对微分方程中的各阶导数进行差分离散，得到差分方程.运用的离散方法有：方法一:用待定系数和泰勒展开进行离散)()()()(d)(d111122iiiiiiiixuxuxuxxu方法二：利用差商逼近导数21122)()(2)()(d)(dhxuxuxuxxuiiiii(2)将(2)带入(1)可以得到)()()()(2)(211uRxfhxuxuxuiiiii,其中)(uRi为无穷小量，这时我们丢弃)(uRi，则有在ix处满足的计算公式：1,...,1)()()(2)(211nixfhxuxuxuiiii，(3)3.根据边界条件，进行边界处理.由(1)可得nuu,0(4)称(3)(4)为逼近(1)的差分方程，并称相应的数值解向量1nU为差分解，iu为)(ixu的近似值.4.最后求解线性代数方程组，得到数值解向量1nU.实验题目：用Lax-Wendroff格式求解方程：.4sin1),1(],1,0[,2sin1)0,(,0),1,0(,02ttuxxxutxuuxt（1）（精确解).2(2sin1txu）数值边值条件分别为：（a）).(20101n0nnnuuhuu（b）.1n0nuu（c）.02-12111n0nnuuu请将计算结果与精确解进行比较。实现算法：1.网格剖分：对求解区域],0()1,0(GT作均匀网格剖分.节点：,jhxjNj,...,1,0,jhtkMk,...,1,0其中空间和时间步长：.,1MTNh2.算法实现将),(1kitxu在节点),(kitx处作泰勒级数展开)(][!2][),(),(32221Otututxutxukikikiki（2）考虑在节点),(kjtx处（1）的微分方程，有：.0xuatu.)()(22222xuatuxaxuattu将上述两式代入（2）式，得)(][2][),(),(322221Oxuaxuatxutxukikikiki对x的一阶、二阶导数用中心差商代替)())],(),([(21][211hOtxutxuhxukikiki)())],(),(2),([(1][211222hOtxutxutxuhxukikikiki代入整理后得到)()()()],(),(2),([2)],(),([2),(),(322211222111OhOhOtxutxutxuhatxutxuhatxutxukikikikikikiki略去误差项，以kiu代替),(kitxu，得到如下差分格式)2(2)(211222111kikikikikikikiuuuhauuhauu（3）（3）式就是Lax-Wendroff格式，其截断误差为)(22hO，节点如图令har||，就得到（1）式的Lax-Wendroff格式的公式)2(2)(2112111kikikikikikikiuuuruuruu（4）（4）式是二阶精度的差分格式.程序代码：function[X,T,U]=advection_fd1d(NS,NT,pde,bd)%WAVE_EQUATION_FD1D利用有限差分方法计算一维双曲线方程%输入参数：%NS整型，空间剖分段数%NT整型，时间剖分段数%pde结构体，带求解的微分方程模型的已知数据，%如边界、初始、系数和右端项等条件.%bd数值边值条件%输出参数:%X长度NS+1的列向量，空间网格剖分%T长度NT+1的行向量，时间网格剖分%U(NS+1)*(NT+1)矩阵，U(:,i)表示第i个时间层网格剖分上的数值解[X,h]=pde.space_grid(NS);[T,tau]=pde.time_grid(NT);N=length(X);M=length(T);U=zeros(N,M);%初值条件U(:,1)=pde.u_initial(X);a=pde.a;r=a*tau/h;%边值条件ifa=0%左边值条件U(1,:)=pde.u_left(T)elseU(end,:)=pde.u_right(T)%右边值条件endfori=2:MU(2:end-1,i)=U(2:end-1,i-1)-r*(U(3:end,i-1)-U(1:end-2,i-1))/2+...r^2*(U(3:end,i-1)-2*U(2:end-1,i-1)+U(1:end-2,i-1))/2;switch(bd)case{'a0'}a0();case{'b'}b();case{'c'}c();otherwisedisp(['Sorry,Idonotknowyour',bd]);endendfunctiona0()U(1,i)=U(1,i-1)-r*(U(2,i-1)-U(1,i-1));endfunctionb()U(1,i)=U(2,i-1);endfunctionc()U(1,i)=2*U(2,i)-U(3,i);endendfunctionpde=model_data()%MODEL_DATA数据模型TI=0;TF=1;SI=0;SF=1;pde=struct('u_exact',@u_exact,'u_initial',@u_initial,...'u_left',@u_left,'u_right',@u_right,'time_grid',...@time_grid,'space_grid',@space_grid,'advection_fd1d_error',@advection_fd1d_error,'a',-2);function[T,tau]=time_grid(NT)T=linspace(TI,TF,NT+1);tau=(TF-TI)/NT;endfunction[X,h]=space_grid(NS)X=linspace(SI,SF,NS+1)'h=(SF-SI)/NS;endfunctionU=u_exact(X,T)[x,t]=meshgrid(X,T);U=1+sin(2*pi*(x+2*t));endfunctionu=u_initial(x)u=1+sin(2*pi*x);endfunctionu=u_right(t)u=1+sin(4*pi*t);endendfunctionshowsolution(X,T,U)%%SHOWSOLUTION以二元函数方式显示数值解%输入参数%X长度为NS+1的列向量，空间网格剖分N%T长度为NT+1的行向量，时间网格剖分M%UN*M矩阵，U(:,i)表示第i个时间层网格部分上的数值解[x,t]=meshgrid(X,T);mesh(x,t,U');xlabel('X');ylabel('T');zlabel('U(X,T)');endfunctionshowvarysolution(X,T,U,UE)%%SHOWVARYSOLUTION显示数值解随着时间的变化%输入参数%X长度为NS+1的列向量，空间网格剖分N%T长度为NT+1的行向量，时间网格剖分M%UN*M矩阵，U(:,i)表示第i个时间层网格部分上的数值解M=size(U,2);figurexlabel('X');ylabel('U');s=[X(1),X(end),min(min(U)),max(max(U))];axis(s);fori=1:Mplot(X,U(:,i));axis(s);pause(0.01);title(['T=',num2str(T(i)),'时刻的温度分布'])End%一维双曲线有限差分方法主测试脚本pde=model_data()[X,T,U]=advection_fd1d(100,200,pde,'a');UE=pde.u_exact(X,T);showvarysolution(X,T,U,UE);%以随时间变化方式显示数值解showsolution(X,T,U);%以二元函数方式显示数值解[X,T,U]=advection_fd1d(100,200,pde,'b');UE=pde.u_exact(X,T);showvarysolution(X,T,U,UE);%以随时间变化方式显示数值解showsolution(X,T,U);%以二元函数方式显示数值解[X,T,U]=advection_fd1d(100,200,pde,'c');UE=pde.u_exact(X,T);showvarysolution(X,T,U,UE);%以随时间变化方式显示数值解showsolution(X,T,U);%以二元函数方式显示数值解