数值分析实验报告之常微分方程数值解汇总
数值分析常微分方程求解实验
实验报告
一、实验目的
解初值问题各种方法比较。
二、实验题目
给定初值问题
⎪⎩⎪⎨⎧=≤<+=,
0)1(,21y x xe x y dx dy x ， 其精确解为)(e e x y x -=，按
（1）改进欧拉法，步长01.0,1.0==h h ；
（2）四阶标准龙格-库塔法，步长1.0=h ；
求在节点)10,...,2,1(1.01=+=k k x k 处的数值解及其误差，比较各个方法的优缺点。
三、实验原理
改进欧拉法程序，四阶标准龙格-库塔法程序。
四、实验内容及结果
五、实验结果分析
实验2中改进欧拉法和四阶标准龙格-库塔法的比较：
结果的第一个ans是x的值与对应的y的值，第二个ans是精确解x的对应值y，第三个ans 是与精确值的误差百分数。
通过误差百分数的比较，可以明显的发现改进欧拉法比四阶标准龙格-库塔法更精确。
常微分方程实验报告
常微分方程实验报告一、实验目的常微分方程是数学分析和实际应用中非常重要的一部分，本次实验的主要目的是通过实际操作和计算，深入理解常微分方程的概念、性质和求解方法，并能够将其应用到实际问题中，提高我们解决数学问题和实际应用问题的能力。
二、实验原理常微分方程是指含有一个自变量和一个未知函数及其导数的等式。
求解常微分方程的方法有很多，常见的有变量分离法、一阶线性方程的求解方法（如常数变易法）、恰当方程的求解方法（通过积分因子）等。
对于一阶常微分方程，形如＼（y' ＋ p(x)y ＝ q(x)＼）的方程，可以使用积分因子＼（e^｛＼int p(x)dx}＼）来求解。
对于可分离变量的方程，形如＼（g(y)dy ＝ f(x)dx\），可以通过分别积分求解。
三、实验内容（一）一阶常微分方程的求解1、求解方程＼（y' ＋ 2xy ＝ 2x\）首先，计算积分因子＼（e^｛＼int 2xdx} ＝ e^｛x^2}＼），然后将方程两边乘以积分因子得到：＼（（ye^｛x^2}）＇＝ 2xe^｛x^2}＼）两边积分可得＼（ye^｛x^2} ＝ e^｛x^2} ＋ C\），解得＼（y ＝1 ＋ Ce^｛x^2}＼）2、求解方程＼（xy' y ＝ x^2\）将方程化为＼（y' ＼frac{y}｛x} ＝ x\），这里＼（p(x) ＝＼frac{1}｛x}＼），积分因子为＼（e^｛＼int ＼frac{1}｛x}dx} ＝＼frac{1}｛x}＼）。
方程两边乘以积分因子得到＼（（＼frac{y}｛x}）＇＝ 1\），积分可得＼（＼frac{y}｛x} ＝ x ＋ C\），即＼（y ＝ x^2 ＋ Cx\）（二）二阶常微分方程的求解1、求解方程＼（y'＇ 2y' ＋ y ＝ 0\）特征方程为＼（r^2 2r ＋ 1 ＝ 0\），解得＼（r ＝ 1\）（二重根），所以通解为＼（y ＝（C_1 ＋ C_2x)e^x\）2、求解方程＼（y'＇＋ 4y ＝ 0\）特征方程为＼（r^2 ＋ 4 ＝ 0\），解得＼（r ＝＼pm 2i\），所以通解为＼（y ＝ C_1\cos(2x) ＋ C_2\sin(2x)＼）（三）应用常微分方程解决实际问题1、考虑一个物体在受到与速度成正比的阻力作用下的运动，其运动方程为＼（m\frac{dv}｛dt} ＝ kv\）（其中＼（m\）为物体质量，＼（k\）为阻力系数），求解速度＼（v\）随时间＼（t\）的变化。
数值分析常微分方程数值解
许多实际问题的数学模型是微分方程或微分方程的定解问题。
如物体运动、电路振荡、化学反映及生物群体的变化等。
常微分方程可分为线性、非线性、高阶方程与方程组等类；线性方程包含于非线性类中，高阶方程可化为一阶方程组。
若方程组中的所有未知量视作一个向量，则方程组可写成向量形式的单个方程。
因此研究一阶微分方程的初值问题⎪⎩⎪⎨⎧=≤≤=0)(),(y a y bx a y x f dxdy， (9-1) 的数值解法具有典型性。
常微分方程的解能用初等函数、特殊函数或它们的级数与积分表达的很少。
用解析方法只能求出线性常系数等特殊类型的方程的解。
对非线性方程来说，解析方法一般是无能为力的，即使某些解具有解析表达式，这个表达式也可能非常复杂而不便计算。
因此研究微分方程的数值解法是非常必要的。
只有保证问题(9-1)的解存在唯一的前提下，研究其数值解法或者说寻求其数值解才有意义。
由常微分方程的理论知，如果(9-1)中的),(y x f 满足条件（1）),(y x f 在区域} ),({+∞<<∞-≤≤=y b x a y x D ，上连续； （2）),(y x f 在上关于满足Lipschitz 条件，即存在常数，使得y y L y x f y x f -≤-),(),(则初值问题(9-1)在区间],[b a 上存在惟一的连续解)(x y y =。
在下面的讨论中，我们总假定方程满足以上两个条件。
所谓数值解法，就是求问题(9-1)的解)(x y y =在若干点b x x x x a N =<<<<= 210处的近似值),,2,1(N n y n =的方法。
),,2,1(N n y n =称为问题(9-1)的数值解，n n x x h -=+1称为由到1+n x 的步长。
今后如无特别说明，我们总假定步长为常量。
建立数值解法，首先要将微分方程离散化，一般采用以下几种方法： (1) 用差商近似导数在问题(9-1)中，若用向前差商hx y x y n n )()(1-+代替)(n x y '，则得)1,,1,0( ))(,()()(1-=≈-+N n x y x f hx y x y n n n n n)(n x y 用其近似值代替，所得结果作为)(1+n x y 的近似值，记为1+n y ，则有 1(,) (0,1,,1)n n n n y y hf x y n N +=+=-这样，问题(9-1)的近似解可通过求解下述问题100(,) (0,1,,1)()n n n n y y hf x y n N y y x +=+=-⎧⎨=⎩(9-2)得到，按式(9-2)由初值经过步迭代，可逐次算出N y y y ,,21。
数值分析常微分方程数值解法
第8页/共105页
➢ 数值积分方法(Euler公式)
设将方程 y=f (x, y)的两端从 xn 到xn+1 求积分, 得
y( xn1) y( xn )
xn1 f ( x, y( x))dx :
xn
xn1 F ( x)dx
xn
用不同的数值积分方法近似上式右端积分, 可以得到计算 y(xn+1)的不同的差分格 式.
h2 2
y''( )
Rn1
:
y( xn1)
yn1
h2 2
y''( )
h2 2
y''( xn ) O(h3 ).
局部截断误差主项
19
第20页/共105页
➢ 向后Euler法的局部截断误差
向后Euler法的计算公式
yn1 yn hf ( xn1, yn1 ), n 0, 1, 2,
定义其局部截断误差为
y 计算 的n递1 推公式，此类计算格式统称为差分格式.
3
第4页/共105页
数值求解一阶常微分方程初值问题
y' f ( x, y), a x b,
y(a)
y0
难点: 如何离散 y ?
➢ 常见离散方法
差商近似导数 数值积分方法 Taylor展开方法
4
第5页/共105页
➢ 差商近似导数(Euler公式)
(0 x 1)
y(0) 1.
解 计算公式为
yn1
yn
hfn
yn
h( yn
2xn ), yn
y0 1.0
n 0, 1, 2,
取步长h=0.1, 计算结果见下表
13
matlab常微分方程的数值解法实验报告
实验四常微分方程的数值解法 指令：[t,y]=ode23(‘fun ’,tspan,yo) 2/3阶龙格库塔方法 [t,y]=ode45(‘fun ’,tspan,yo) 4/5阶龙格库塔方法 [t,y]=ode113(‘fun ’,tspan,yo) 高阶微分方程数值方法其中fun 是定义函数的文件名。
该函数fun 必须以为dx 输出量，以t,y 为输入量。
tspan=[t0 tfina]表示积分的起始值和终止值。
yo 是初始状态列向量。
考虑到初始条件有00d , (0)0,d d , (0)0.d SSI S S tI SI I I I tββμ⎧=-=>⎪⎪⎨⎪=-=≥⎪⎩ (5.24) 这就是Kermack 与McKendrick 的SIR 仓室模型. 方程(5.24)无法求出()S t 和()I t 的解析解.我们先做数值计算。
Matlab 代码为：function dy=rigid(t,y) dy=zeros(2,1); a=1; b=0.3;dy(1)=a*y(1).*y(2)-b*y(1); dy(2)=-a*y(1).*y(2);ts=0:.5:50; x0=[0.02,0.98];[T,Y]=ode45('rigid',ts,x0); %plot(T,Y(:,1),'-',T,Y(:,2),'*') plot(Y(:,2),Y(:,1),'b--') xlabel('s') ylabel('i')任务：1 画出i （t ），2分析各参数的影响例57：求解两点边值问题：0)5(,0)1(,32==='-''y y x y y x 。
（注意：相应的数值解法比较复杂）。
y=dsolve('x*D2y-3*Dy=x^2','y(1)=0,y(5)=0','x') ↙ y =-1/3*x^3+125/468+31/468*x^4例：用数值积分的方法求解下列微分方程 π21''2t y y -=+设初始时间t0=0；终止时间tf=3*pi ；初始条件0|',0|00====x x y y 。
微分方程数值解法实验报告
微分方程数值解法实验报告班级：姓名：学号：日期：一、实验目的1、熟悉微分方程（组）数值解的Euler算法，改进的Euler算法和Runge-Kutta算法，利用matlab软件实现微分方程数值解法来求解具体试题；2、虽然求解常微分方程有各种各样的解析解，但解析方法只能用来求解一些特殊类型的方程，通常它们无法求出解析解，而需要数值方法来近似求解。
因此产生了常微分方程初值问题的数值计算方法，常微分方程数值解法是通过计算机便捷的求解近似值。
二、基本理论及背景1、在数值求解常微分方程中，主要有有限差分计算和有限元计算两大类方法，其中在有限差分计算方法中有一类方法称为龙格－库塔(Runge_Kutta)方法。
四阶的龙格－库塔方法为最佳的计算格式。
2、参考三中的代码，分别用Euler算法，改进的Euler算法和Runge-Kutta 算法实现微分方程（组）的数值求解，完成下列题目：三、算法设计及实现1、算法设计，通过Euler算法，改进的Euler算法和Runge-Kutta三种算法来实现微分方程（组）的数值求解；2、程序文件及功能清单：(1) Euler Method:function [x,y]=EulerDSolve(f,ab,y0,h)x=(ab(1):h:ab(2))';n=length(y0);y=zeros(length(x),n);y(1,:)=y0';for k=2:length(x)y(k,:)=y(k-1,:)+h*feval(f,x(k-1),y(k-1,:)')';end；(2) Improved Euler Method:function [x,y]=MEulerDSolve(f,ab,y0,h)x=(ab(1):h:ab(2))';n=length(x);y=zeros(n,length(y0));y(1,:)=y0';for k=2:nyp=y(k-1,:)+h*feval(f,x(k-1),y(k-1,:));yc=y(k-1,:)+h*feval(f,x(k),yp);y(k,:)=(yp+yc)/2;end(3) Runge-Kutta Method:function [TOut,YOut]=Runge_Kutta(f,ab,y0,h)TOut=(ab(1):h:ab(2))';n=length(TOut);YOut=zeros(n,length(y0));YOut(1,:)=y0';for k=2:nx=TOut(k-1); y=YOut(k-1,:)';K1=feval(f,x,y);K2=feval(f,x+h/2,y+K1*h/2);K3=feval(f,x+h/2,y+K2*h/2);K4=feval(f,x+h,y+K3*h);YOut(k,:)=(y+(K1+2*K2+2*K3+K4)*h/6)';end四、实验步骤1、打开MATLAB软件，新建 *.m文件，在m文件的窗口中编辑Euler算法的函数程序，另建一m文件，编辑自己改进的Euler算法的函数程序，再新建一m文件，在窗口中编辑Runge-Kutta算法的函数程序，并全部保存在指定的文件夹下；2、将MATLAB软件的工作页面的工具栏下的目标文件指向指定的文件夹；3、分别调用上述三种算法的函数，实现微分方程（组）的数值求解完成给定的实验题目；4、输出结果和初步分析说明（见附页）。
【清华】实验4-常微分方程数值解
实验4 常微分方程的数值解【实验目的】1. 掌握用MA TLAB 软件求微分方程初值问题数值解的方法；2. 通过实例学习用微分方程模型解决简化的实际问题；3. 了解欧拉方法和龙格-库塔方法的基本思想和计算公式，及稳定性等概念。
【实验内容】 1.题目2.小型火箭初始重量为1400kg ，其中包括1080kg 燃料。
火箭竖直向上发射时燃料燃烧率为18kg/s ，由此产生32000N 的推力，火箭引擎在燃烧用尽时关闭。
设火箭上升时空气阻力正比于速度的平方，比例系数为0.4kg/m ，求引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点时的高度和加速度，并画出高度、速度、加速度随时间变化的图形。
模型及其求解 假设在上升过程中，重力加速度不随高度的变化，恒为g =9.8m/s 2。
（1）分析引擎关闭前，火箭质量为m ，高度为h ，速度为v ，加速度为a ，阻力为f ，则由题意可知：()140018m m t t ==-，dhv dt=，20.4f v = 由牛顿第二定律可得：2320000.49.8140018F dv F f mg v a dt m m t---====--总由此得到一个常微分方程组：dhv dt= 2320000.49.8140018dv v dt t -=-- 初值为v =0m/s,h =0m ；条件为1080kg60s 18kg/st ≤=根据常微分方程组的初值问题，在MA TLAB 中计算数值解，记(1)x h =，(2)x v=T ((1),(2))x x x =。
源程序为：由MA TLAB绘图得到：由MA TLAB计算得到从开始上升到关闭引擎瞬间的情况如下表：表1.所以，引擎关闭瞬间,火箭的高度为h=12190m，速度v=267.26m/s，加速度a=0.91701m/s^2。
（2）引擎关闭后，20.49.8320dv vadt==--，dhvdt=；初值为（1）中的末态数值，条件是60st>,且速度v>0，用循环语句找出符合上述条件的时间范围[60,T]。
常微分方程的数值解法
常微分方程的数值解法1. 引言常微分方程是自变量只有一个的微分方程，广泛应用于自然科学、工程技术和社会科学等领域。
由于常微分方程的解析解不易得到或难以求得，数值解法成为解决常微分方程问题的重要手段之一。
本文将介绍几种常用的常微分方程的数值解法。
2. 欧拉方法欧拉方法是最简单的一种数值解法，其具体步骤如下：- 将自变量的区间等分为n个子区间；- 在每个子区间上假设解函数为线性函数，即通过给定的初始条件在每个子区间上构造切线；- 使用切线的斜率（即导数）逼近每个子区间上的解函数，并将其作为下一个子区间的初始条件；- 重复上述过程直至达到所需的精度。
3. 改进的欧拉方法改进的欧拉方法是对欧拉方法的一种改进，主要思想是利用两个切线的斜率的平均值来逼近每个子区间上的解函数。
具体步骤如下： - 将自变量的区间等分为n个子区间；- 在每个子区间上构造两个切线，分别通过给定的初始条件和通过欧拉方法得到的下一个初始条件；- 取两个切线的斜率的平均值，将其作为该子区间上解函数的斜率，并计算下一个子区间的初始条件；- 重复上述过程直至达到所需的精度。
4. 二阶龙格-库塔方法二阶龙格-库塔方法是一种更为精确的数值解法，其基本思想是通过近似计算解函数在每个子区间上的平均斜率。
具体步骤如下： - 将自变量的区间等分为n个子区间；- 在每个子区间上计算解函数的斜率，并以该斜率的平均值近似表示该子区间上解函数的斜率；- 利用该斜率近似值计算下一个子区间的初始条件，并进一步逼近解函数；- 重复上述过程直至达到所需的精度。
5. 龙格-库塔法（四阶）龙格-库塔法是目前常用的数值解法之一，其精度较高。
四阶龙格-库塔法是其中较为常用的一种，其具体步骤如下：- 将自变量的区间等分为n个子区间；- 在每个子区间上进行多次迭代计算，得到该子区间上解函数的近似值；- 利用近似值计算每个子区间上的斜率，并以其加权平均值逼近解函数的斜率；- 计算下一个子区间的初始条件，并进一步逼近解函数；- 重复上述过程直至达到所需的精度。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
数学与计算科学学院 实 验 报 告
实验项目名称 常微分方程数值解 所属课程名称 数值方法B 实 验 类 型 验证 实 验 日 期 2013.11.11
班 级 学 号 姓 名 成 绩 1
一、实验概述： 【实验目的】
1.掌握求解常微分方程的欧拉法； 2.掌握求解常微分方程的预估校正法； 3.掌握求解常微分方程的经典的四阶龙格库塔法； 4.能用C语言或MATLAB将上述三种算法用程序运行出来； 5.将算法实例化，并得出三种算法的相关关系，如收敛性、精度等； 6.附带书中例题的源程序见附录1。
【实验原理】 1．欧拉格式 （1）显式欧拉格式：1(,)nnnnyyhfxy
局部截断误差：22211()()()()22nnnhhyxyyyxoh （2）隐式欧拉格式：111(,)nnnnyyhfxy 局部截断误差：2211()()()2nnnhyxyyxoh
2．预估校正法 预估：1(,)nnnnyyhfxy 校正：111[(,)(,)]2nnnnnnhyyfxyfxy 2
统一格式：1[(,)(,(,))]2nnnnnnnnhyyfxyfxhyhfxy
平均化格式：11(,),(,),1().2pnnncnnpnpcyyhfxyyyhfxyyyy
3．四阶龙格库塔方法的格式（经典格式）
112341213243(22),6(,),(,),22(,),22(,).nnnnnnnnnnhyyKKKKKfxyhhKfxyKhhKfxyKKfxhyhK
【实验环境】 1.硬件环境： HP Microsoft 76481-640-8834005-23929 HP Corporation Intel(R) Core(TM) I5-2400 CPU @ 3.10GHz 3.09GHz,3.16GB的内存
2.软件环境： Microsoft Windows XP 3
Professional 版本 2002 Service Pack 3
二、实验内容： 【实验方案】
方案一： 用欧拉法，预估校正法，经典的四阶龙格库塔方法求解下列ODE问题：
例题：在区间【0，1】上以h=0.1用欧拉法，预估校正法,经典的四阶龙格库塔
法求解微分方程 dy/dx=-y+x+1，初值y（0）=1；其精确解为y=x+exp(-x),且将计算结果与精确解进行比较，对三个算法的收敛性的进行分析比较。
方案二： 用欧拉法，预估校正法, 经典的四阶龙格库塔方法求解初值问题 dy/dx=xxey，初值y（0）=1； 将计算结果与精确解为21(2)2xxe
比较在区
间[0,1]上分别取步长h=0.1; 0.05时进行计算。对三个算法的收敛性进行分析比较，
【实验过程】（实验步骤、记录、数据、分析） 注：以下图形是通过Excel表格处理数据得出，并未通过MATLAB编程序所得！ 1、1(0)1dyyxdxy 4
由题可知精确解为：xyxe，当x=0时，y(x)=0。
h=0.1 表1 h=0.1时三个方法与精确值的真值表
图1 h=0.1时三个方法走势图 h=0.05（此时将源程序中i的范围进行扩大，即for(i=0;i<20;i++)）
步长 Euler法 预估校正法 经典四阶库 精确值 0.1 1.010000 1.005000 1.004838 1.249080 0.2 1.029000 1.019025 1.018731 1.055455 0.3 1.056100 1.041218 1.040818 1.091217 0.4 1.090490 1.070802 1.070320 1.131803 0.5 1.131441 1.107076 1.106531 1.176851 0.6 1.178297 1.149404 1.148812 1.226025 0.7 1.230467 1.197211 1.196586 1.279016 0.8 1.287421 1.249975 1.249329 1.335536 0.9 1.348678 1.307228 1.306570 1.395322 1.0 1.413811 1.368541 1.367880 1.458127 5
表2 h=0.05时三个方法与精确值的真值表 步长 Euler法 预估校正法 经典四阶库 精确值 0.05 1.002500 1.001250 1.001229 1.011721 0.10 1.007375 1.004877 1.004837 1.024908 0.15 1.014506 1.010764 1.010708 1.039504 0.20 1.023781 1.018802 1.018731 1.055455 0.25 1.035092 1.028885 1.028801 1.072710 0.30 1.048337 1.040915 1.040818 1.091217 0.35 1.063421 1.054795 1.054688 1.110931 0.40 1.080250 1.070436 1.070320 1.131801 0.45 1.098737 1.087752 1.087628 1.153791 0.50 1.118800 1.106662 1.106531 1.176851 0.55 1.140360 1.127087 1.126950 1.200942 0.60 1.163342 1.148954 1.148812 1.226025 0.65 1.187675 1.172193 1.172046 1.252062 0.70 1.213291 1.196736 1.196585 1.279016 0.75 1.240127 1.222520 1.222367 1.306852 0.80 1.268121 1.249485 1.249329 1.335536 0.85 1.297215 1.277572 1.277415 1.365037 0.90 1.327354 1.306728 1.306570 1.395322 0.95 1.358486 1.336900 1.336741 1.426362 1.00 1.390562 1.368039 1.367880 1.458127
图2 h=0.05时三个方法走势图 6
2、(0)1xdyxeydxy 由题可知精确解为：21(2)2xxe，当x=0时，y(x)=0。 h=0.1 表3 h=0.1时三个方法与精确值的真值表 步长 Euler法 预估校正法 经典四阶库 精确值 0.1 0.900000 0.909625 0.909428 0.929533 0.2 0.819249 0.835927 0.835593 0.872564 0.3 0.754433 0.776081 0.775655 0.826822 0.4 0.702726 0.727671 0.727189 0.790348 0.5 0.661726 0.688636 0.688127 0.761457 0.6 0.629396 0.657225 0.656711 0.738709 0.7 0.604018 0.631957 0.631453 0.720874 0.8 0.584147 0.611582 0.611100 0.706908 0.9 0.568575 0.595050 0.594599 0.695927 1.0 0.556297 0.581487 0.581072 0.687191 7
图3 h=0.1时三个方法走势图 h=0.05（此时将源程序中i的范围进行扩大，即for(i=0;i<20;i++)） 表4 h=0.05时三个方法与精确值的真值表 步长 Euler法 预估校正法 经典四阶库 精确值 0.05 0.950000 0.952452 0.952427 0.962924 0.10 0.904904 0.909474 0.909428 0.929533 0.15 0.864284 0.870670 0.870606 0.899511 0.20 0.827741 0.835671 0.835592 0.872564 0.25 0.794908 0.804137 0.804047 0.848419 0.30 0.765447 0.775755 0.775655 0.826822 0.35 0.739043 0.750232 0.750125 0.807538 0.40 0.715407 0.727302 0.727189 0.790348 0.45 0.694272 0.706715 0.706599 0.775050 0.50 0.675394 0.688245 0.688126 0.761457 0.55 0.658546 0.671682 0.671561 0.749397 0.60 0.643519 0.656830 0.656710 0.738709 0.65 0.630124 0.643514 0.643395 0.729247 0.70 0.618185 0.631570 0.631453 0.720874