微分方程数值解法课程试验题目
《微分方程数值解》课程论文题目1
《微分方程数值解》课程论文题目1、给出Dirichlet 边值问题：()()()(),u x q x u x f x a x b''-+=<< (),()u a u b αβ==和两点边值问题(())()()(),,d d u d u L u p x r x q x u f x a x b d x d x d x =-++=<< (),()u a u b αβ==在均匀网格下的差分格式，编写程序，给出误差阶及误差图。
2、考虑守恒型微分方程（I ）：(())()(),d du Lu p x q x u f x a x b dx dx =-+=<< (),()u a u b αβ==给出其直接差分格式和积分插值法差分格式，采用积分插值法时，数值积分分别利用中矩形公式和梯形公式，编写程序，给出误差阶及误差图。
3、考虑守恒型微分方程（II ）：(())()(),d du Lu p x q x u f x a x b dx dx =-+=<< 01()(),u a u a αα'=+01()().u b u b ββ'=+给出其直接差分格式和积分插值法差分格式，对于边界条件采取直接微分法和积分插值法，编写程序，给出误差阶及误差图。
4、考虑Poisson 方程：(,),(,)u f x y x y G -∆=∈ |(,),u x y αΓ=其中G 是xy 平面上一有界区域，其边界Γ为光滑曲线。
给出方程的五点差分格式和九点差分格式，给出其截断误差，编写程序，给出误差阶及误差图。
5、考虑Poisson 方程：(,),(,)u f x y x y G -∆=∈ |(,),u x y αΓ=其中G 是xy 平面上一圆域、环形域或扇形域，其边界Γ为光滑曲线。
给出方程极坐标形式的差分格式，编写程序，给出误差阶及误差图。
6、考虑Laplace 方程：0,(,),u x y G -∆=∈ |(,),u x y αΓ=其中G 是xy 平面上一有界区域，其边界Γ为光滑曲线。
微分方程数值解作业1
微分方程数值解法---实验报告专业：信息与计算科学班级：计算071姓名：梁成保学号：3070811027西安理工大学2010-6-12第一章 常微分方程初值问题数值解法实验一1、实验题目试用(a)欧拉格式(b)中点格式 (c)预报—校正格式 (d)经典四级四阶R-K 格式 编程计算方程：(]()⎪⎩⎪⎨⎧=∈+=112,1222y x y x x y dx dy 2、程序#include<iostream.h> #include<stdlib.h> #include<math.h> const int N=11;double fund(double x,double y); void Euler(double a,double h,double y[]); void Center(double a,double h,double y[]); void YuJiao(double a,double h,double y[]); void SjSj(double a,double h,double y[]); void Adams(double a,double h,double y[]); void main() { double a,b,h,y[N]; int i; char option;cout<<"请输入定义域区间[a,b]：\n";cin>>a>>b;cout<<"请输入初值y[0]：\n";cin>>y[0];h=(b-a)/(N-1);cout<<"a:欧拉法\nb:中心法\nc:预报-校正格式\n";cout<<"d:经典四级四阶R-K格式\n";cout<<"e:Adams预报-修正格式\nf:退出\n";label: cout<<"请选择算法：";cin>>option;switch(option){case 'a': Euler(a,h,y);break;case 'b': Center(a,h,y);break;case 'c': YuJiao(a,h,y);break;case 'd': YuJiao(a,h,y);break;case 'e': Adams(a,h,y);break;case 'f': exit(1);break;default : {cout<<"选择有错,重新选择！\n";goto label;}}cout<<"计算结果为：\n xn\t\tyn"<<endl;for(i=0;i<N;i++)cout<<" "<<a+i*h<<"\t\t"<<y[i]<<endl;}void Euler(double a,double h,double y[]){int i;for(i=1;i<N;i++){y[i]=y[i-1]+h*fund(a+(i-1)*h,y[i-1]);}}void Center(double a,double h,double y[]){int i;double w;for(i=1;i<N;i++){w=fund(a+(i-1)*h,y[i-1]);y[i]=y[i-1]+h*fund(a+(i-1)*h/2,w*h/2+y[i-1]);}}void YuJiao(double a,double h,double y[]){int i;double sun;double y_Begin[N]={y[0]};for(i=1;i<N;i++){y_Begin[i]=y_Begin[i-1]+h*fund(a+(i-1)*h,y_Begin[i-1]);}for(i=1;i<N;i++){while(fabs(y_Begin[i]-sun)>0.0001){sun=y[i-1]+h/2.0*(fund(a+(i-1)*h,y[i-1])+fund(a+i*h,y_Begin[i]));y_Begin[i]=sun;}y[i]=sun;}}void SjSj(double a,double h,double y[]){int i;double k1,k2,k3,k4;for(i=1;i<N;i++){k1=fund(a+(i-1)*h,y[i-1]);k2=fund(a+(i-1)*h+h/2,y[i-1]+h*k1/2);k3=fund(a+(i-1)*h+h/2,y[i-1]+h*k2/2);k4=fund(a+(i-1)*h+h,y[i-1]+h*k3);y[i]=y[i-1]+h/6*(k1+2*k2+2*k3+k4);}}void Adams(double a,double h,double y[]){int i;double me,x[N];double k1,k2,k3,k4;for(i=0;i<N;i++){x[i]=a+i*h;y[i]=y[0];}for(i=1;i<4;i++){k1=fund(a+(i-1)*h,y[i-1]);k2=fund(a+(i-1)*h+h/2,y[i-1]+h*k1/2);k3=fund(a+(i-1)*h+h/2,y[i-1]+h*k2/2);k4=fund(a+(i-1)*h+h,y[i-1]+h*k3);y[i]=y[i-1]+h/6*(k1+2*k2+2*k3+k4);}for(i=4;i<N;i++){me=y[i-1]+h/24*(55*fund(x[i-1],y[i-1])-59*fund(x[i-2],y[i-2])+37*fund(x[i-3],y[i-3])-9*fund(x[i-4],y[i-4]));while(fabs(me-y[i])>0.0001){y[i]=y[i-1]+h/24*(9*fund(x[i],me)+19*fund(x[i-1],y[i-1])-5*fund(x[i-2],y[i-2])+fund(x[i-3],y[i-3])); me=y[i]; }}}double fund(double x,double y) { double s;s=y/(2*x)+x*x/2/y; return s; }3、试验结果及分析(i)运算结果(ii)结果分析在使用欧拉法数值求解过程中，其计算过程非常简单，即由0y 可直接计算出1y ，由1y可直接计算出,,2 y 无需用迭代方法求解任何方程，就可求出近似解n y 。
第10章实验九常微分方程初值问题数值解
第10章实验九常微分方程初值问题数值解第10章实验九常微分方程初值问题数值解实验目的：掌握Euler 方法、改进的Euler 方法、Runge-Kutta 方法，并能在MATLAB 中应用这些方法来数值计算10.1 Euler 方法微分方程在科学和工程中常用于建立数学模型，众所周知，求解微分方程的解析解常常较困难，这时就需要数值近似解。
对于如下的初值问题：],[)(),(0b a x y a y y x f y ∈==' （10.1）我们并非要找一个满足初值问题的可微函数，而是生成点集{}k k y x ,(，并以这些点作为近似值即k k y x y ≈)(。
如何构造近似值k y 呢？首先选择点的横坐标，一般将区间[a,b]划为n 个子区间，选择横坐标点为：n a b h n k kh a x k -==+=其中,,,1,0, 值h 称为步长，然后近似求解为：1,1,0),(),(10001-=+=+=+n k y x hf y y y x hf y y k k k k （10.2）这种求近似解的方法被称为Euler 方法。
例10.1 对于微分方程 1/3(1)1y ty y '=??=?，用Euler 方法，分别取h=0.1,0.05,0.01,...求解，计算到(2),y 将数值解与精确解[]2/323/)2()(+=t t y 比较。
解：（1）写出欧拉显示方法源程序代码(,)f t y 的.m 文件：function E=euler(f,a,b,ya,M)h=(b-a)/M;t=zeros(1,M+1);y=zeros(1,M+1);yy=zeros(1,M+1);z=zeros(1,M+1);t=a:h:b; % 选好步长y(1)=ya;for j=1:My(j+1)=y(j)+h*feval(f,t(j),y(j)); % 计算近似值yk endyy=((t.^2+2)/3).^(3/2); % 原方程的解析解z=yy-y; % 解析解与近似解之差fprintf('\n Extended Trapezoidal Rule\n');fprintf('\n t y yy Error \n');E=[t',y',z'];plot(t,y,t,yy,'r',t,z)（2）写出函数(,)f t y的.m文件：function f=w1_(t,y)f=t*y^(1/3);（3）写出执行函数命令：E=euler(f,a,b,ya,M)其中f表示求导函数，a,b分别为左右端点，ya为初始值，M为步幅数，并将a,b,ya,M分别取定为1，2，1，M（10，20，100），就得到如下执行语句E=euler2('w1_',1,2,1,10) ;E=euler2('w1_',1,2,1,20) ;E=euler2('w1_',1,2,1,100)。
数学实验课程设计常微分方程数值解
数学实验报告1.题目：某容器盛满水后，底端直径为d0的小孔开启（如图1），根据水力学知识，当水面高度为h时，谁从小孔中流出的速度为v=0.6*（g*h）^0.5（其中g为重力加速度，0.6问哦小孔收缩系数）1）若容器为倒圆锥形（如图1），现测得容器高和上底直径都为1.2m，小孔直径d为3cm，为水从小孔中流完需要多少时间；2min时水面高度是多少。
2）若容器为倒葫芦形（如图2），现测得容器高1.2m，小孔直径d为3cm，由底端（记x=0）向上每隔0.1m测出容器的直径D（m）如表1所示，问水从小孔中流完需要多少时间；2min时水面高度是多少。
X/m00.10.20.30.40.50.60.70.80.9 1.0 1.1 1.2D/m0.030.050.080.140.190.330.450.680.981.11.21.131.02.分析：由题知，水从小孔中流出，不仅与容器有关，还与水流速度v=0.6*（2*g*h）^0.5有关。
第一小题容器是圆锥形，比较规则，但是由于水不断从小孔流出，容器中水的高度是不断变化的，水流速度没有一定的公式，所以要用到微积分解决。
由（1）知，水面直径等于水深。
水深为h时，流量为0.6(π/4)d^2*(gh)^0.5，0.6*(g*h)^(0.5)*π*(d0/2)^2*dt=π/4*h^2*dh则水深下降dh 所需时间 ：dt=-[(π/4)h^2*dh]/[0.6(π/4)d^2*(gh)^0.5]=-[h^1.5*dh]/[0.6d^2*(g)^0.5]水深由1.2m 至0定积分得水从小孔流完的时间：T （其中已知d=0.03m ，g=9.8m*s(-2）对于第二问：设两分钟（120S ）后水深为X m ，由dt=-[(π/4)h^2*dh]/[0.6*(π/4)*d^2*(gh)^0.5]=-[h^1.5*dh]/[0.6d^2*(g)^0.5]则263.93-120 =X^2.5/[1.5*d^2*(g)^0.5]以d=0.03m ，g=9.8m*s(-2代入上式得 水深：X第二小题容器为倒葫芦形，比较不规则，比较复杂，不仅要考虑水不断从小孔流出，容器中水的高度是不断变化的，水流速度没有一定的公式，所以要用到微积分解决，还要注意表1的倒葫芦形的不断变化，水深的高度变化是不规则的但仍可以用微积分。
微分方程数值解习题课
微分方程初值问题数值解习题课一、应用向前欧拉法和改进欧拉法求由如下积分2xt y e dt -=⎰所确定的函数y 在点x =0.5,1.0,1.5的近似值。
解：该积分问题等价于常微分方程初值问题2'(0)0x y e y -⎧=⎪⎨=⎪⎩其中h=0.5。
其向前欧拉格式为2()100ih i i y y he y -+⎧=+⎪⎨=⎪⎩改进欧拉格式为22()2(1)10()20ih i h i i h y y ee y --++⎧=++⎪⎨⎪=⎩将两种计算格式所得结果列于下表二、应用4阶4步阿达姆斯显格式求解初值问题'1(0)1y x y y =-+⎧⎨=⎩00.6x ≤≤ 取步长h=0.1.解：4步显式法必须有4个起步值，0y 已知，其他3个123,,y y y 用4阶龙格库塔方法求出。
本题的信息有：步长h=0.1；结点0.1(0,1,,6)i x ih i i === ；0(,)1,(0)1f x y x y y y =-+==经典的4阶龙格库塔公式为11234(22)6i i hy y k k k k +=++++1(,)1i i i i k f x y x y ==-+121(,)0.05 1.0522i i i i hk hk f x y x y k =++=--+232(,)0.05 1.0522i i i i hk hk f x y x y k =++=--+433(,)0.1 1.1i i i i k f x h y hk x y k =++=--+算得1 1.0048375y =,2 1.0187309y =,3 1.0408184y =4阶4步阿达姆斯显格式1123(5559379)24i i i i i i hy y f f f f +---=+-+-1231(18.5 5.9 3.70.90.24 3.24)24i i i i i y y y y y i ---=+-+++ 由此算出4561.0703231, 1.1065356, 1.1488186y y y ===三、用Euler 方法求()'1,0101x y e y x x y =-++≤≤=问步长h 应该如何选取，才能保证算法的稳定性？解：本题(),1xf x y e y x =-++ (),0,01x y f x y e x λ'==-<≤≤ 本题的绝对稳定域为111x h he λ+=-<得02x he <<，故步长应满足02,00.736he h <<<<四、 求梯形方法111[(,)(,)]2k k k k k k hy y f x y f x y +++=++的绝对稳定域。
实验 4 常微分方程的数值解 几道解方程的题——重要
G、空气阻力 f 和推力 F。其中推力 F 为常量：F 32000N；空气阻力正比于速度
的平方，
， 为比例系数，由模型假设，k 认为是常量，即 k=0.4kg/m；重
力 G （M‐mt）g，M 1400kg 为火箭初始重量，m 18kg/s 为燃料燃烧率，t
为火箭发射后的飞行时间，g 为重力加速度，由模型假设，认为 g=9.8N/kg 为常 量。 开始发射后火箭具有向上的加速度，速度不断增大，阻力也随着不断增大； 由牛顿第二定律 可知：
题目 2
【题目 1】 （课本习题第四章第 2 题） 小型火箭初始重量为 1400kg，其中包括 1080kg 燃料。火箭竖直向上发射时 燃料燃烧率为 18kg/s，由此产生 32000N 的推力，火箭引擎在燃料用尽时关闭。 设火箭上升时空气阻力正比于速度的平方，比例系数为 0.4kg/m，求引擎关闭瞬 间火箭的高度、速度、加速度，及火箭到达最高点时的高度和加速度，并画出高 度、速度、加速度随时间变化的图形。 【模型假设】 基于对题目的分析及求解需要，可作以下假设： 1.假设火箭在飞行过程中重力加速度的数值不变，恒为 9.8m/s2（由于火箭
30
20
10
0
-10
-20
-30
-40
0
20
40
60
80
100
120
但如果步长选得不一样(横坐标都从
微分方程数值解实习题
function y=f2(t,u)
y=t^3*u^3-t*u;
%Adams二步外插+Euler方法，例2
n=input('n=');%n等分
T=input('T=');%时间上界
H=input('H=');%第一个小区间H等分
h=T/n;%步长
hh=h/H;%小步长
u=zeros(n+1,1);
Columns 31 through 36
0.1865 0.1910 0.1955 0.2000 0.2043 0.2086
Columns 37 through 42
0.2128 0.2170 0.2211 0.2252 0.2292 0.2331
Columns 43 through 48
0.2370 0.2409 0.2447 0.2484 0.2521 0.2558
for m=1:n-1
u1(m+2)=u1(m+1)+(1/12)*h*(5*f1((m+2)*h,u1(m+2))+8*f1((m+1)*h,u1(m+1))-f1(m*h,u1(m)));
end
uu=zeros(n+1,1);
uu(1)=1;
for m=2:n+1
uu(m)=1/3*exp(m*h)+2/3*exp(-2*m*h);
Columns 57 through 63
0.2872 0.2906 0.2939 0.2972 0.3005 0.3037 0.3070
Columns 64 through 70
0.3102 0.3134 0.3166 0.3198 0.3230 0.3261 0.3293
实验5常微分方程的数值解
实验5 常微分方程的数值解概要：将装满放射性废物的圆桶扔到水深300ft 的海底，圆桶体积55gal ，装满废料的桶重为527.436lbf ，在海中浮力为470.327lbf 。
此外，下沉时受到的阻力与速度成正比，比例系数为0.08lbf/s 。
实验发现当圆桶速度超过40ft/s 时，就会因与海底冲撞而破裂。
要求：（1）建立解决上述问题的微分方程模型（2）用数值和解析两种方法求解微分方程，并回答谁赢得了官司。
模型建立由牛顿第二定律可列出圆桶下沉速度的微分方程：dv G F f G F bv dt m m ----==其中G 为圆桶重量，F 为浮力，b 为下沉阻力与速度关系的比例系数。
换算到国际单位制，dept=300*0.3048=91.4400 海深（m ）ve=40*0.3048=12.1920 速度极限（超过就会使圆筒碰撞破裂）(m/s) G=527.436*0.4536*9.8＝2344.6 圆筒重量（N ） F=470.327*0.4536*9.8＝2090.7 浮力(N)m=527.436*0.4536＝239.24 圆筒质量（kg ） b=0.08*0.4536*9.8/0.3048＝1.1667 比例系数(Ns/m) 模型求解 一．求数值解Matlab 程序如下： sd.m:function dx=sd(t,x,G,F,m,b) dx=[(G-F-b*x)/m];%微分方程sddraw.m: clear;G=527.436*0.4536*9.8;%圆筒重量（N ） F=470.327*0.4536*9.8;%浮力(N)m=527.436*0.4536;%圆筒质量（kg ）b=0.08*0.4536*9.8/0.3048%比例系数(Ns/m) h=0.1;%所取时间点间隔ts=[0:h:2000];%粗略估计到时间2000 x0=0;%初始条件opt=odeset('reltol',1e-3,'abstol',1e-6);%相对误差1e-6，绝对误差1e-9 [t,x]=ode45(@sd,ts,x0,opt,G ,F,m,b);%使用5级4阶龙格—库塔公式计算 %[t,x]%输出t,x(t),y(t)plot(t,x,'-'),grid%输出v(t)的图形 xlabel('t'); ylabel('v(t)');%用辛普森公式对速度积分求出下沉深度 T=20;%估计20s 以内降到海底for i=0:2:10*T%作图时间间隔为0.2 y=x(1:(i+1)); k=length(y);a1=[y(2:2:k-1)];s1=sum(a1); a2=[y(3:2:k-1)];s2=sum(a2);z4((i+2)/2)=(y(1)+y(k)+4*s1+2*s2)*h/3;%辛普森公式求深度 endi=[0:2:10*T]; figure;de=300.*0.3048.*ones(5*T+1,1);%海深ve=40.*0.3048*[1 1];%速度极限值（超过就会使圆筒碰撞破裂）plot(x(i+1),z4',x(i+1),de,ve,[0 z4(5*T+1)]);%作出速度－深度图线，同时画出海深和速度要求grid;gtext('dept'),gtext('Vmax');xlabel('v');ylabel('dept(v)');figure;plot(i/10,z4');%作出时间－下降深度曲线grid;xlabel('t');ylabel('dept(t)');求解结果如下图：速度—时间曲线：可以看到经过足够长的时间后，如若桶没有落到海底，它的速度会趋于常值。
常微分方程数值解练习
h¯ ∈ C | h¯j/j!| < 1
j=0
4
1、用待定系数法确定系数，使下面方法之阶尽可能高，并求局部截断
误差的主项：
(1)、un+1 = α1un + α2un−1 + hβ0fn+1.
(2)、un+1 = un + h(β1fn + β2fn−1).
(3)、un+1 = α2un1 + h(β0fn+1 + β1fn + β2fn−1) 2、讨论下列方法的相容性，稳定性和绝对稳定性：
(1)、un+1 = −4un + 5un−1 + h(4fn + 2fn−1),
(2)、un+1
=
1 8
(9un
− un−2) +
3 8
h(fn+1
+
2fn
− fn−1)
3、求α，使线性多步法
un+2
−
(1
+
α)un+1
+
αun
=
h 2
[(3
−
α)fn+1
−
(1
+
α)fn]
是相容的和稳定的。
4、证明三阶Adams公式
2
1、研究求解：
du dt
=
f (t, u), 0
<
t
≤
T,
u(0) = u0
的方法
un+1 = un + h[θf (tn, un) + (1 − θ)f (tn+1, un+1)]
微分方程数值解课程设计题目
课程设计题目1.计算并画出二级二阶、三级三阶和四级四阶Runge-Kutta 方法的绝对稳定区域。
2.用古典隐式格式并结合追赶法计算抛物型方程2201,00.01(,0)sin()01(0,)(1,)00u u x t t x u x x x u t u t t π⎧∂∂= <<<<⎪∂∂⎪⎪= ≤≤⎨⎪== ≥⎪⎪⎩的近似解。
其中0.1,0.001h k ==。
并与解析解2sin t u ex ππ−=比较。
3.用五点差分格式,1,1,,1,1,21(4)0l m l m l m l m l m l m U U U U U U h +−+−◊=+++−= 求解Dirichlet 问题 {}2222220(,),(,)0,1ln (1)u u x y x y x y xy u x y ∂Ω⎧∂∂+= ∈ΩΩ=<<⎪∂∂⎨⎪⎡⎤=++⎣⎦⎩ 其中，1.0=Δ=Δy x ，分别用Jacobi 迭代和Gauss －Seidel 迭代求解，误差为310−，分析两方法的收敛速度，将计算结果画出图形。
4.方程2004,012,1(,0)1,1u u x t t x x u x x ∂∂⎧+= <<<<⎪∂∂⎪⎨≤⎧⎪=⎨⎪>⎩⎩其中0.1h =，0.01k =。
分别用Lax-Friedrichs 格式、C-I-R 格式和Lax-Wendroff 。
格式求解上述微分方程初边值问题。
并画出图形比较计算结果。
5.对初边值问题2222000101,011sin ,018|0,010t t x x u u x t tx u x x u x tu u t π====⎧∂∂= <<<<⎪∂∂⎪⎪=≤≤⎪⎨⎪∂=≤≤⎪∂⎪⎪== 0≤≤1⎩ 利用显式差分格式：)2(211211nm n m n m n m n m n m U U U r U U U −+−++−=+− 求近似解。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
计算实验课
微分方程数值解法数值计算实验题目
一、常微分方程部分：
1．使用四阶Runge-Kutta 方法求解如下初值问题的近似解，并将结果与实际值进行比较。
2．使用四阶Adams 预估校正算法（PECP 和PMECME 方案），初始值用四阶Runge-Kutta 方法提供，并将结果与实际值进行比较。
()2
1u t u -+='，32≤≤t ，()12=u ；精度510-=ε，5.0=h 。
实际解11u t t
=+-。
t
u
u +='1，21≤≤t ，()21=u ；精度510-=ε，2.0=h 。
实际解2ln +=t t u 。
二、偏微分方程部分：
1．用有限差分法求解如下Poisson 方程
(),cos3sin ,u x y x y π-∆=，0x π<<，10<<y
边界条件为： ()(),0,10,u x u x ==01x ≤≤； ()()0,,0,x x u y u y π==10≤≤y 取1,h N
π
=
和21
,h N
=
作矩形剖分，网格节点为1i x ih =，2j y jh =，i ,j =0,1,…,N 。
差分格式为
1,,1,,1,,11222cos3sin i j i j i j i j i j i j i j u u u u u u x y h h π+-+--+-+⎛⎫-+= ⎪⎝⎭
,1,2,,1i j N =-L 边界条件为： 00,0,,,i iN u u i N ===L
01,1,,1,j j u u j N ==-L 1,1,,1,N j N j u u j N -==-L 结果与精确解()()
1
2,9cos3sin u x y x y π
π-=+进行比较。
求解方案：依次令4,8,16,32N =，取6位小数计算。
用消元法求解，并就(),,44j i j x y π⎛⎫
= ⎪⎝⎭
，,1,2,3i j =处列出差分解与精确解。
其次，就N =32，0.25，0.5，0.75及i =0，2，4，…，30，32画出差分解曲
线。
2．用向前、向后或Crank-Nicolson 算法求解一维抛物型方程的初边值问题：
22
sin ,u u
t t x ∂∂=+∂∂，01x <<，t <0 ()0,x u t =()1,0x u t =，t <0 (),0cos u x x π=，01x << 精确解为：()()2
,cos 1cos t u x t e x t ππ-=+-，设空间步长1
h J
=，时间步长0τ>，t k =k t ，网比2
r h
τ
=。
（一） 向前差分格式的计算方案：
111
2
2sin k k k k k
j j
j j j k u u u u u t h
τ
++---+=
+ 1,2,,1j J =-L ,1,2,,k N =L
边值条件为 j =0，
100
101
2
2sin k k k k k
k u u u u u t h
τ
+---+=+，11k n u u -=; j =J ，1112
2sin k k k k k
J J
J J J k u u u u u t h
τ
++---+=+，11k n J J u u -+=; 初值条件为
(),0cos ,
0,1,2,,j j u x x j J π==L
a) 取1,40h =
1,3200τ= 此时12r =，计算到时间层32001t =； b) 取1,80h = 1,12800τ= 此时1
2r =，计算到时间层128001t =；
c) 取1,80h = 1
,3200
τ= 此时2r =，观察计算结果；
（二） 向后差分格式的计算方案：
111
12
2sin k k k k k
j j
j j j k u u u u u t h τ
++-+--+=
+ 1,2,,1j J =-L ,1,2,,k N =L
边值条件为 j =0，
1111
00
101
12
2sin k k k k k k u u u u u t h
τ
++++-+--+=+，1111k k u u ++-=; j =J ，1112
2sin k k J J
J J J k u u u u u t h
τ
++---+=
+，11
11k k J J u u ++-+=; 初值条件为
(),0cos ,
0,1,2,,j j u x x j J π==L
a) 取1,40h =
1,1600τ= 此时1r =，计算到时间层16001t =； b) 取1,80h = 1
,3200
τ= 此时2r =，计算到时间层32001t =；
（三） 六点对称差分格式的计算方案：
()()()(
)
()1111111112
1
1
2sin sin 22
k k j j
k k k k k k j j j j j j k k u u u u u u u u t t h τ
++++++--+-=
+-++++
+ 1,2,,1j J =-L ,1,2,,k N =L 边值条件为
j =0，J 列Crank-Nicolson 格式，其中
111111k k k k u u u u ++--+=+; 111111k k k k J J J J u u u u ++----+=+
初值条件为
(),0cos ,
0,1,2,,j j u x x j J π==L
a) 取1,40h =
1
,1600τ= 此时1r =，计算到时间层16001t =； b) 取1,80h = 1,3200
τ= 此时2r =，计算到时间层32001t =； 将以上三种数值方法的结果与精确解列表作比较，其中,1,,44
j j
x j ==L 。
二维抛物型方程的初边值问题*：
用六点对称差分格式，ADI 法，预校法和LOD 法求解二维抛物型方程的初边值问题：
2
22224,u u u t x y -⎛⎫∂∂∂=+ ⎪∂∂∂⎝⎭
()()(),0,10,1x y G ∈=⨯，t <0 ()0,,u y t =()1,,0u y t =，01y <<，t <0 (),0,y u x t =(),1,0y u x t =，01x <<，t <0 (),,0cos sin u x y y x ππ=， 精确解为：()2
8
,,sin cos t
u x y t e
x y πππ-
=。
设x j =jh (j =0,1,…,J ), y k =kh (k =0,1,…,K ), t n =n t (n =0,1,…,N )，差分解为n jk u ，则边值条件为
00n n
k Jk u u ==，k =0,1,…,K ；
01n n j j u u =，1n n jK jK u u -=，j =0,1,…,J
取空间步长12140h h h ===
，时间步长11600τ=，网比21r h
τ
==，用六点对称差分格式，ADI 法，预校法和LOD 法分别计算到时间层t =1。
3．微分方程数值解法书第181页的实习题1、2。