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

《偏微分方程数值和matlab实验报告(8)》由会员上传分享，免费在线阅读，更多相关内容在工程资料-天天文库。

1、WORD格式.可编辑专业知识.整理分享WORD格式.可编辑实验题目：用Lax-Wendroff格式求解方程：（1）（精确解）数值边值条件分别为：（a）（b）（c）请将计算结果与精确解进行比较。实现算法：1.网格剖分：对求解区域作均匀网格剖分.节点：其中空间和时间步长：2.算法实现将在节点处作泰勒级数展开（2）考虑在节点处（1）的微分方程，有：将上述两式代入（2）式，得对的一阶、二阶导数用中心差商代替专业知识.整理分享WORD格式.可编辑代入整理后得到略去误差项，以代替，得到如下差分格式（3）（3）式就是Lax-Wendroff格式，其截断误差为，节

2、点如图令，就得到（1）式的Lax-Wendroff格式的公式（4）（4）式是二阶精度的差分格式.程序代码：function[X,T,U]=advection_fd1d(NS,NT,pde,bd)%WAVE_EQUATION_FD1D利用有限差分方法计算一维双曲线方程专业知识.整理分享WORD格式.可编辑%输入参数：%NS整型，空间剖分段数%NT整型，时间剖分段数%pde结构体，带求解的微分方程模型的已知数据，%如边界、初始、系数和右端项等条件.%bd数值边值条件%输出参数:%X长度NS+1的列向量，空间网格剖分%T长度NT+1的行向量，时间网格剖分

3、%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

4、-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()专业知识.整理分享WORD格式.可编辑U(1,i)=U(1,i-1)-r*(U(2,i-1)-U(1,i-1));endfunctio

5、nb()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_f

6、d1d_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*

7、pi*x);endfunctionu=u_right(t)u=1+sin(4*pi*t);endend专业知识.整理分享WORD格式.可编辑functionshowsolution(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)');en

8、dfunctionshowvarysolution(X,T,U,UE)%%SHOWVARYSOLUTION显示数值解随着时