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

实现算法： 1. 网格剖分：对求解区域 G (0,1) (0,T ] 作均匀网格剖分 .节点：x jjh , j 0,1,..., Nt k jh ,k0,1,..., M其中空间和时间步长：h1 , T .NM2. 算法实现将 u(x i , t k 1 ) 在节点 ( x i ,t k ) 处作泰勒级数展开u(x i , t k 1) u(x i , t k )[u k 22u kO ( 3)]i[t 2 ]it2!考虑在节点 ( x j , t k ) 处（ 1）的微分方程，有：u a u 0.t x2u( a u a ( u22ut 2) ) ax 2.txx t将上述两式代入（ 2）式，得u( x i ,t k 1 ) u( x i ,t k ) a [ u] i k a 2 2 [ 2 u 2 ]i k O ( 3 )x 2 x对 x 的一阶、二阶导数用中心差商代替（ 1）（ 2）u k12[x ]i2h[( u( x i1, t k )u( x i 1, t k ))]O(h)[2u k12u( x i ,t k )u( x i 1, t k ))]2) x2 ]i h2 [( u(x i 1 ,t k )O (h代入整理后得到u( x i ,t k 1 ) u( x i ,t k )a[ u( x i 1 ,t k ) u( x i 1 ,t k )] 2ha 22O( h2 ) 2 h2 )O( 3 ) 2[ u( x i 1 ,t k )2u( x i , t k )u( x i1 , t k )]O (2h略去误差项，以u i k代替 u( x i ,t k ) ，得到如下差分格式a(u i k122u i k 1u i k u i k1 )a2 (u i k12u i k u i k1)（ 3）2h2h（ 3）式就是 Lax-Wendroff 格式，其截断误差为O(2h2 ) ，节点如图| a |，就得到（1）式的 Lax-Wendroff格式的公式令 r hu i k 1u i k r(u i k1u i k1 )r 2(u i k12u i k u i k1)（ 4）22（ 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_left',@u_left, 'u_right',@u_right,,@u_initial,'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);%以二元函数方式显示数值解'c' );[X,T,U]=advection_fd1d(100,200,pde,UE=(X,T);showvarysolution(X,T,U,UE);%以随时间变化方式显示数值解showsolution(X,T,U);%以二元函数方式显示数值解。