./實(shí)驗(yàn)題目：用Lax-Wendroff格式求解方程：〔1〔精確解數(shù)值邊值條件分別為：請(qǐng)將計(jì)算結(jié)果與精確解進(jìn)行比較。實(shí)現(xiàn)算法：網(wǎng)格剖分：對(duì)求解區(qū)域作均勻網(wǎng)格剖分.節(jié)點(diǎn)：其中空間和時(shí)間步長(zhǎng)：算法實(shí)現(xiàn)將在節(jié)點(diǎn)處作泰勒級(jí)數(shù)展開〔2考慮在節(jié)點(diǎn)處〔1的微分方程,有：將上述兩式代入〔2式,得對(duì)的一階、二階導(dǎo)數(shù)用中心差商代替代入整理后得到略去誤差項(xiàng),以代替,得到如下差分格式〔3〔3式就是Lax-Wendroff格式,其截?cái)嗾`差為,節(jié)點(diǎn)如圖令,就得到〔1式的Lax-Wendroff格式的公式〔4〔4式是二階精度的差分格式.程序代碼：function[X,T,U]=advection_fd1d<NS,NT,pde,bd>%WAVE_EQUATION_FD1D利用有限差分方法計(jì)算一維雙曲線方程%輸入?yún)?shù)：%NS整型,空間剖分段數(shù)%NT整型,時(shí)間剖分段數(shù)%pde結(jié)構(gòu)體,帶求解的微分方程模型的已知數(shù)據(jù),%如邊界、初始、系數(shù)和右端項(xiàng)等條件.%bd數(shù)值邊值條件%輸出參數(shù):%X長(zhǎng)度NS+1的列向量,空間網(wǎng)格剖分%T長(zhǎng)度NT+1的行向量,時(shí)間網(wǎng)格剖分%U<NS+1>*<NT+1>矩陣,U<:,i>表示第i個(gè)時(shí)間層網(wǎng)格剖分上的數(shù)值解[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數(shù)據(jù)模型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以二元函數(shù)方式顯示數(shù)值解%輸入?yún)?shù)%X長(zhǎng)度為NS+1的列向量,空間網(wǎng)格剖分N%T長(zhǎng)度為NT+1的行向量,時(shí)間網(wǎng)格剖分M%UN*M矩陣,U<:,i>表示第i個(gè)時(shí)間層網(wǎng)格部分上的數(shù)值解[x,t]=meshgrid<X,T>;mesh<x,t,U'>;xlabel<'X'>;ylabel<'T'>;zlabel<'U<X,T>'>;endfunctionshowvarysolution<X,T,U,UE>%%SHOWVARYSOLUTION顯示數(shù)值解隨著時(shí)間的變化%輸入?yún)?shù)%X長(zhǎng)度為NS+1的列向量,空間網(wǎng)格剖分N%T長(zhǎng)度為NT+1的行向量,時(shí)間網(wǎng)格剖分M%UN*M矩陣,U<:,i>表示第i個(gè)時(shí)間層網(wǎng)格部分上的數(shù)值解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>>,'時(shí)刻的溫度分布']>End%一維雙曲線有限差分方法主測(cè)試腳本pde=model_data<>[X,T,U]=advection_fd1d<100,200,pde,'a'>;UE=pde.u_exact<X,T>;showvarysolution<X,T,U,UE>;%以隨時(shí)間變化方式顯示數(shù)值解showsolution<X,T,U>;%以二元函數(shù)方式顯示數(shù)值解[X,T,U]=advection_fd1d<100,200,pde,'b'>;UE=pde.u_exact<X,T>;showvarysolution<X,T,U,UE>;%以隨時(shí)間變化方式顯示數(shù)值解showsolution