偏微分方程數(shù)值及matlab實驗報告(共8頁)

1、精選優(yōu)質(zhì)文檔-傾情為你奉上實驗題目：用Lax-Wendroff格式求解方程： （1） （精確解）數(shù)值邊值條件分別為：（a）（b）（c）請將計算結果與精確解進行比較。實現(xiàn)算法：1. 網(wǎng)格剖分：對求解區(qū)域作均勻網(wǎng)格剖分.節(jié)點： 其中空間和時間步長：2. 算法實現(xiàn)將在節(jié)點處作泰勒級數(shù)展開 （2）考慮在節(jié)點處（1）的微分方程，有：將上述兩式代入（2）式，得對的一階、二階導數(shù)用中心差商代替代入整理后得到略去誤差項，以代替，得到如下差分格式 （3）（3）式就是Lax-Wendroff格式，其截斷誤差為，節(jié)點如圖 令，就得到（1）式的Lax-Wendroff格式的公式 （4）（4）式是二階精度的差分格式.程

2、序代碼：function X,T,U = advection_fd1d (NS ,NT ,pde,bd)% WAVE_EQUATION_FD1D 利用有限差分方法計算一維雙曲線方程% 輸入?yún)?shù)：% NS 整型，空間剖分段數(shù)% NT 整型，時間剖分段數(shù)% pde 結構體，帶求解的微分方程模型的已知數(shù)據(jù)，% 如邊界、初始、系數(shù)和右端項等條件.% bd 數(shù)值邊值條件% 輸出參數(shù):% X 長度 NS+1 的列向量，空間網(wǎng)格剖分% T 長度 NT+1 的行向量，時間網(wǎng)格剖分% U (NS+1)*(NT+1) 矩陣，U(:,i) 表示第 i 個時間層網(wǎng)格剖分上的數(shù)值解 X,h = pde.space_gr

3、id(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;% 邊值條件if a >= 0 % 左邊值條件 U(1,:) = pde.u_left(T)else U(end,:) = pde.u_right(T) %右邊值條件end for i = 2:M U(2:end -1,i) =U(2:end-1,i-1)-r*(U(3:end,i-1)-U(1:end-2,i-1)/2+.

4、 r2*(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(); otherwise disp('Sorry, I do not know your ', bd); end end function a0() U(1,i)=U(1,i-1)-r*(U(2,i-1)-U(1,i-1); end function b() U(1,i)=U(2,i-1); end function c(

5、) U(1,i)=2*U(2,i)-U(3,i); end endfunction pde = 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,'

6、advection_fd1d_error',advection_fd1d_error,'a',-2); function T,tau = time_grid(NT) T = linspace(TI,TF,NT+1); tau = (TF-TI)/NT; end function X,h = space_grid(NS) X = linspace(SI,SF,NS+1)' h = (SF-SI)/NS; end function U = u_exact(X,T) x,t=meshgrid(X,T); U = 1+sin(2*pi*(x+2*t); end func

7、tion u = u_initial (x) u = 1+sin(2*pi*x); end function u = u_right(t) u =1+sin(4*pi*t); endend function showsolution (X,T,U)% SHOWSOLUTION 以二元函數(shù)方式顯示數(shù)值解% 輸入?yún)?shù)% X 長度為 NS +1 的列向量，空間網(wǎng)格剖分N% T 長度為 NT +1 的行向量，時間網(wǎng)格剖分M% U N*M 矩陣，U(:,i)表示第i個時間層網(wǎng)格部分上的數(shù)值解x,t = meshgrid (X,T);mesh (x,t,U');xlabel ('X'

8、;);ylabel ('T');zlabel ('U(X,T)');endfunction showvarysolution (X,T,U,UE)% SHOWVARYSOLUTION 顯示數(shù)值解隨著時間的變化% 輸入?yún)?shù)% X 長度為 NS +1 的列向量，空間網(wǎng)格剖分N% T 長度為 NT +1 的行向量，時間網(wǎng)格剖分M% U N*M 矩陣，U(:,i)表示第i個時間層網(wǎng)格部分上的數(shù)值解M = size (U ,2) ;figurexlabel ('X');ylabel ('U');s = X(1) ,X(end),min(mi

9、n(U),max(max(U);axis (s);for i = 1:M plot (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);%以隨時間變化方式顯示數(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ù)值解showsolution(X,T,U);%以二元函數(shù)方式顯示數(shù)值解