微分方程数值解法课程试验题目

《微分方程数值解》课程论文题目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 平面上一有界区域，其边界Γ为光滑曲线。

微分方程数值解法---实验报告专业：信息与计算科学班级：计算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章实验九常微分方程初值问题数值解实验目的：掌握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 +++=++的绝对稳定域。

实验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)');求解结果如下图：速度—时间曲线：可以看到经过足够长的时间后，如若桶没有落到海底，它的速度会趋于常值。

课程设计题目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 −+−++−=+− 求近似解。