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

./实验题目：用Lax-Wendroff格式求解方程：〔1〔精确解数值边值条件分别为：请将计算结果与精确解进行比较。实现算法：网格剖分：对求解区域作均匀网格剖分.节点：其中空间和时间步长：算法实现将在节点处作泰勒级数展开〔2考虑在节点处〔1的微分方程,有：将上述两式代入〔2式,得对的一阶、二阶导数用中心差商代替代入整理后得到略去误差项,以代替,得到如下差分格式〔3〔3式就是Lax-Wendroff格式,其截断误差为,节点如图令,就得到〔1式的Lax-Wendroff格式的公式〔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