偏微分方程數值和matlab實驗報告

./實驗題目：用Lax-Wendroff格式求解方程：〔1〔精確解數值邊值條件分別為：請將計算結果與精確解進行比較。實現算法：網格剖分：對求解區(qū)域作均勻網格剖分.節(jié)點：其中空間和時間步長：算法實現將在節(jié)點處作泰勒級數展開〔2考慮在節(jié)點處〔1的微分方程,有：將上述兩式代入〔2式,得對的一階、二階導數用中心差商代替代入整理后得到略去誤差項,以代替,得到如下差分格式〔3〔3式就是Lax-Wendroff格式,其截斷誤差為,節(jié)點如圖令,就得到〔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