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

实验题目：用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; % 边值条件 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+... 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(); 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() U(1,i)=2*U(2,i)-U(3,i); end end function pde = 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; end function [X,h] = space_grid(NS)