偏微分方程数值及matlab实验报告(共8页).doc

《偏微分方程数值及matlab实验报告(共8页).doc》由会员分享，可在线阅读，更多相关《偏微分方程数值及matlab实验报告(共8页).doc（8页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、精选优质文档-倾情为你奉上实验题目：用Lax-Wendroff格式求解方程： （1） （精确解）数值边值条件分别为：（a）（b）（c）请将计算结果与精确解进行比较。实现算法：1. 网格剖分：对求解区域作均匀网格剖分.节点： 其中空间和时间步长：2. 算法实现将在节点处作泰勒级数展开 （2）考虑在节点处（1）的微分方程，有：将上述两式代入（2）式，得对的一阶、二阶导数用中心差商代替代入整理后得到略去误差项，以代替，得到如下差分格式 （3）（3）式就是Lax-Wendroff格式，其截断误差为，节点如图 令，就得到（1）式的Lax-Wendroff格式的公式 （4）（4）式是二阶精度的差分格式.程

2、序代码：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_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+. r2*

4、(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 endfunction pd

5、e = 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

6、 = (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 function 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)% SHO

7、WSOLUTION 以二元函数方式显示数值解% 输入参数% X 长度为 NS +1 的列向量，空间网格剖分N% T 长度为 NT +1 的行向量，时间网格剖分M% U N*M 矩阵，U(:,i)表示第i个时间层网格部分上的数值解x,t = meshgrid (X,T);mesh (x,t,U);xlabel (X);ylabel (T);zlabel (U(X,T);endfunction showvarysolution (X,T,U,UE)% SHOWVARYSOLUTION 显示数值解随着时间的变化% 输入参数% X 长度为 NS +1 的列向量，空间网格剖分N% T 长度为 NT +1

8、 的行向量，时间网格剖分M% U N*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);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);U

9、E=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(X,T,U);%以二元函数方式显示数值解 X,T,U=advection_fd1d(100,200,pde,c);UE=pde.u_exact(X,T);showvarysolution(X,T,U,UE);%以随时间变化方式显示数值解showsolution(X,T,U);%以二元函数方式显示数值解专心-专注-专业