偏微分方程数值实验报告八实验题目：利用有限差分法求解.0)1(,0)1(),()()(==-=+''-u u x f x u x u 真解为)1()(22x e x u x -=-实现算法：对于两点边值问题,)(,)(,,d 22βα==∈=-b u a u l x f dxu(1)其中),(b a l =f b a ),(<为],[b a l =上的连续函数，βα，为给定常数.其相应的有限差分法的算法如下：1.对求解区域做网格剖分，得到计算网格.在这里我们对区间l 均匀剖分n 段，每个剖分单元的剖分步长记为nab h -=.2.对微分方程中的各阶导数进行差分离散，得到差分方程.运用的离散方法有：方法一:用待定系数和泰勒展开进行离散)()()()(d )(d 111122++--++≈i i i i i i i i x u x u x u x x u ααα方法二：利用差商逼近导数21122)()(2)()(d )(d h x u x u x u x x u i i i i i -++-≈(2)将(2)带入(1)可以得到)()()()(2)(211u R x f h x u x u x u i i i i i +=+---+,其中)(u R i 为无穷小量，这时我们丢弃)(u R i ，则有在i x 处满足的计算公式：1,...,1)()()(2)(211-==+---+n i x f hx u x u x u i i i i ，(3)3.根据边界条件，进行边界处理.由(1)可得βα==n u u ,0(4)称(3)(4)为逼近(1)的差分方程，并称相应的数值解向量1-n U 为差分解，i u 为)(i x u 的近似值.4.最后求解线性代数方程组，得到数值解向量1-n U .实验题目：用Lax-Wendroff 格式求解方程：.4sin 1),1(],1,0[,2sin 1)0,(,0),1,0(,02t t u x x x u t x u u x t ππ+=∈+=>∈=- （1） （精确解).2(2sin 1t x u ++=π） 数值边值条件分别为： （a ）).(20101n 0nn nu u hu u -+=+τ （b ）.1n 0nu u =（c ）.02-12111n 0=++++n n u u u请将计算结果与精确解进行比较。

实现算法： 1.网格剖分：对求解区域],0()1,0(G T ⨯=作均匀网格剖分. 节点： ,jh x j = N j ,...,1,0= ,jh t k = M k ,...,1,0= 其中空间和时间步长：.,1MT N h ==τ 2.算法实现将),(1+k i t x u 在节点),(k i t x 处作泰勒级数展开)(][!2][),(),(32221τττO tu t u t x u t x u ki k i k i k i +∂∂+∂∂+=+ （2）考虑在节点),(k j t x 处（1）的微分方程，有：.0=∂∂-∂∂xu a t u .)()(22222x u a t u x a x u a t t u ∂∂=∂∂∂∂-=∂∂-∂∂=∂∂ 将上述两式代入（2）式，得)(][2][),(),(322221τττO xu a x u a t x u t x u ki k i k i k i +∂∂+∂∂-=+对x 的一阶、二阶导数用中心差商代替)())],(),([(21][211h O t x u t x u hx u k i k i k i +-=∂∂-+ )())],(),(2),([(1][211222h O t x u t x u t x u hx u k i k i k i k i ++-=∂∂-+ 代入整理后得到)()()()],(),(2),([2)],(),([2),(),(322211222111τττττO h O h O t x u t x u t x u ha t x u t x u ha t x u t x u k i k i k i k i k i k i k i ++++-+--=-+-++略去误差项，以ki u 代替),(k i t x u ，得到如下差分格式)2(2)(211222111k i k i k i ki k i k ik iu u u ha u u h a u u---+++-+--=ττ （3） （3）式就是Lax-Wendroff 格式，其截断误差为)(22h O +τ，节点如图 令ha r τ||=，就得到（1）式的Lax-Wendroff 格式的公式 )2(2)(2112111k i k i k i ki k i k ik iu u u r u u r u u---+++-+--= （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] = (NS);[T,tau] = (NT);N = length(X); M = length(T);U = zeros(N,M);% 初值条件U(:,1) = (X);a = ;r = a*tau/h;% 边值条件if a >= 0 % 左边值条件U(1,:) = (T)elseU(end,:) = (T) %右边值条件endfor i = 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, I do not know your ', bd]);endendfunction a0()U(1,i)=U(1,i-1)-r*(U(2,i-1)-U(1,i-1));endfunction b()U(1,i)=U(2,i-1);endfunction c()U(1,i)=2*U(2,i)-U(3,i);endendfunction 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;endfunction [X,h] = space_grid(NS)X = linspace(SI,SF,NS+1)'h = (SF-SI)/NS;endfunction U = u_exact(X,T)[x,t]=meshgrid(X,T);U = 1+sin(2*pi*(x+2*t));endfunction u = u_initial (x)u = 1+sin(2*pi*x);endfunction u = u_right(t)u =1+sin(4*pi*t);endendfunction showsolution (X,T,U)%% SHOWSOLUTION 以二元函数方式显示数值解% 输入参数% 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 的行向量，时间网格剖分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:Mplot (X,U(:,i));axis (s);pause ;title ([ 'T=',num2str(T(i)),'时刻的温度分布'])End% 一维双曲线有限差分方法主测试脚本pde=model_data()[X,T,U]=advection_fd1d(100,200,pde,'a');UE=(X,T);showvarysolution(X,T,U,UE);%以随时间变化方式显示数值解showsolution(X,T,U);%以二元函数方式显示数值解[X,T,U]=advection_fd1d(100,200,pde,'b');UE=(X,T);showvarysolution(X,T,U,UE);%以随时间变化方式显示数值解showsolution(X,T,U);%以二元函数方式显示数值解[X,T,U]=advection_fd1d(100,200,pde,'c');UE=(X,T);showvarysolution(X,T,U,UE);%以随时间变化方式显示数值解showsolution(X,T,U);%以二元函数方式显示数值解。