常微分方程的数值解法实验报告
常微分方程实验报告
常微分方程实验报告一、实验目的常微分方程是数学分析和实际应用中非常重要的一部分，本次实验的主要目的是通过实际操作和计算，深入理解常微分方程的概念、性质和求解方法，并能够将其应用到实际问题中，提高我们解决数学问题和实际应用问题的能力。
二、实验原理常微分方程是指含有一个自变量和一个未知函数及其导数的等式。
求解常微分方程的方法有很多，常见的有变量分离法、一阶线性方程的求解方法（如常数变易法）、恰当方程的求解方法（通过积分因子）等。
对于一阶常微分方程，形如＼（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\）的变化。
常微分方程实验报告
常微分方程实验报告《常微分方程》综合性实验实验报告实验班级05应数（3）学生姓名江晓荣学生学号200530770314指导老师方平华南农业大学理学院应用数学系实验微分方程在数学建模中的应用及数值解的求法一、实验目的1．了解常微分方程的基本概念。
2．常微分方程的解了解析解和数值解。
3．学习、掌握MA TLAB 软件有关求解常微分方程的解析解和数值解的有关命令。
4. 掌握微分方程在数学建模中的应用。
二、实验内容1．用MA TLAB 函数dsolve 符号求解常微分方程的通解和特解。
2．用MA TLAB 软件数值求解常微分方程。
三、实验准备1．用MA TLAB 求常微分方程的解析解的命令用MA TLAB 函数dsolve 求常微分方程()(,,,,,,)0n F x y y y y y ''''''= （7.1）的通解的主要调用格式如下：S=dsolve('eqn', 'var')其中输入的量eqn 是改用符号方程表示的常微分方程(,,,2,)0F x y Dy D y Dny = ，导数用D 表示，2阶导数用D2表示，以此类推。
var 表示自变量，默认的自变量为t 。
输出量S 是常微分方程的解析通解。
如果给定常微分方程（7.1）的初始条件()00010(),(),,()n n y x a y x a y x a '=== ，则求方程（7.1）的特解的主要调用格式如下：S=dsolve('eqn', 'condition1 ',…'conditonn ','var')其中输入量eqn ，var 的含义如上，condition1,…conditonn 是初始条件。
输出量S 是常微分方程的特解。
2．常微分方程的数值解法除常系数线性微分方程可用特征根法求解、少数特殊方程可用初等积分法求解外，大部分微分方程无解析解，应用中主要依靠数值解法。
微分方程数值解法实验报告
微分方程数值解法实验报告班级：姓名：学号：日期：一、实验目的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、输出结果和初步分析说明（见附页）。
微分方程数值解实验报告
微分方程数值解实验报告实验目的：掌握微分方程数值解的基本方法，能够利用计算机编程求解微分方程。
实验原理：微分方程是自然科学与工程技术中常见的数学模型，它描述了变量之间的关系及其随时间、空间的变化规律。
解微分方程是研究和应用微分方程的基础，但有很多微分方程无法找到解析解，只能通过数值方法进行求解。
本实验采用欧拉方法和改进的欧拉方法求解微分方程的初值问题：$$\begin{cases}\frac{dy}{dt}=f(t,y)\\y(t_0)=y_0\end{cases}$$其中，$f(t,y)$是给定的函数，$y(t_0)=y_0$是已知的初值条件。
欧拉方法是最基本的数值解法，其步骤如下：1.给定$t_0$和$y_0$2.计算$t_{i+1}=t_i+h$，其中$h$是步长3. 计算$y_{i+1}=y_i+hf(t_i,y_i)$4.重复步骤2、3直到达到终止条件改进的欧拉方法是对欧拉方法进行改进，通过利用函数$y(t)$在$t+\frac{1}{2}h$处的斜率来更准确地估计$y_{i+1}$，其步骤如下：1.给定$t_0$和$y_0$2.计算$t_{i+1}=t_i+h$，其中$h$是步长3. 计算$y_*=y_i+\frac{1}{2}hf(t_i,y_i)$4. 计算$y_{i+1}=y_i+hf(t_i+\frac{1}{2}h,y_*)$5.重复步骤2、3、4直到达到终止条件实验步骤：1.编写程序实现欧拉方法和改进的欧拉方法2.给定微分方程和初值条件3.设置步长和终止条件4.利用欧拉方法和改进的欧拉方法求解微分方程5.比较不同步长下的数值解与解析解的误差6.绘制误差-步长曲线，分析数值解的精度和收敛性实验结果：以一阶常微分方程$y'=3ty+t$为例，给定初值$y(0)=1$，取步长$h=0.1$进行数值求解。
利用欧拉方法求解微分方程得到的数值解如下：\begin{array}{cccc}t & y_{\text{exact}} & y_{\text{Euler}} & \text{误差} \\ \hline0.0&1.000&1.000&0.000\\0.1&1.035&1.030&0.005\\0.2&1.104&1.108&0.004\\0.3&1.212&1.217&0.005\\0.4&1.360&1.364&0.004\\0.5&1.554&1.559&0.005\\0.6&1.805&1.810&0.005\\0.7&2.131&2.136&0.005\\0.8&2.554&2.560&0.006\\0.9&3.102&3.107&0.006\\1.0&3.807&3.812&0.005\\\end{array}利用改进的欧拉方法求解微分方程得到的数值解如下：\begin{array}{cccc}t & y_{\text{exact}} & y_{\text{Improved Euler}} & \text{误差} \\\hline0.0&1.000&1.000&0.000\\0.1&1.035&1.035&0.000\\0.2&1.104&1.103&0.001\\0.3&1.212&1.211&0.001\\0.4&1.360&1.358&0.002\\0.5&1.554&1.552&0.002\\0.6&1.805&1.802&0.003\\0.7&2.131&2.126&0.005\\0.8&2.554&2.545&0.009\\0.9&3.102&3.086&0.015\\1.0&3.807&3.774&0.032\\\end{array}误差-步长曲线如下：实验结论：通过对比欧拉方法和改进的欧拉方法的数值解与解析解的误差，可以发现改进的欧拉方法具有更高的精度和收敛性。
微分方程数值方法实验报告
微分方程数值方法实验报告微分方程数值方法实验报告一、实验目的1、了解Euler法及梯形法的基本原理，能用其解决常微分方程初值问题，并把计算结果用图形表示出来。
2、理解4阶RK法基本计算步骤，能编程实现算法并解决相关常微分方程初值问题。
3、了解MATLAB主要功能和基本特征，熟悉MATLAB操作环境。
掌握MATLAB常用函数的使用以及图形处理。
二、实验题目对于初值问题u’=u,u(0)=1,在区间[0,1]上，用Euler法，梯形法及RK方法进行计算，取步长h=0.1,0.2,0.5，试比较(1)用同样步长，三种方法中哪一个精度最好;(2)对同一种方法一不同步长进行计算，哪一个结果最好。
三、实验内容1、步长为h=0.1时，用三种方法计算题目1)、MATLAB程序Euler法:>> a=0;b=1;h=0.1;N=(b-a)/h;y=zeros(N+1,1);y(1)=1;x=a:h:b;>> for i=2:N+1y(i)=y(i-1)+h*y(i-1);end求得:y = (1 1.1 1.21 1.331 1.4641 1.6105 1.7716 1.9487 2.1436 2.3579 2.5937)’梯形法:>> z=zeros(N+1,1);>> z(1)=1;>> for i=2:N+1v1=z(i-1);t=z(i-1)+h*v1;v2=t;z(i)=z(i-1)+h/2*(v1+v2);end1求得:z = (1 1.105 1.2205 1.3476 1.4873 1.641 1.8101 1.99622.2008 2.4258 2.6734)’RK法:>> w=zeros(N+1,1);>> w(1)=1;>> for i=2:N+1K1=w(i-1);K2=w(i-1)+K1*h/2;K3=w(i-1)+K2*h/2;K4=w(i-1)+K3*h;w(i)=w(i-1)+h*(K1+2*K2+2*K3+K4)/6;end求得:w =(1 1.1052 1.2214 1.3499 1.4918 1.6487 1.8221 2.01382.2255 2.4596 2.7183)’>> plot(x,y)>> hold on>> plot(x,z,':')>> plot(x,w,'--')>> plot(x,exp(x),'*')>> title('相同步长下三种方法与准确解的对比')>> legend('Euler法','梯形法','四阶RK法','准确解') 2)、图形对比相同步长下三种方法与准确解的对比2.82.62.42.221.81.6Euler法1.4梯形法四阶RK法1.2准确解1 00.10.20.30.40.50.60.70.80.9123)、结果分析由图得出三种方法中四阶RK法经确度最高，梯形法次之，Euler法精确度最差。
常微分方程上机实验报告
实验名称
常微分方程初值问题的数值解法
实验时间
2012年3月13日
姓名
XXX
班级
XXX
学号
XXX
成绩
一、实验目的及内容
1．MATLAB基本操作：矩阵相乘、矩阵LU分解及矩阵分解为另一个矩阵和其转置相乘、解矛盾方程、符号求解微分方程。
2．MATLAB数据调用及绘图。
二、相关背景知识介绍
（4）符号求解微分方程，求 。
2.MATLAB数据调用及绘图：
三．代码
1.(1)A=[3 4;5 6]
B=[0 8;1 6]
A*B
(2)A=[2 3;4 5]
[L U]=lu(A)
B=Chol(A)
(3)A=[2 8;3 4]
b=[6;2]
X=A\b
(4)dsolve('Dy=x','x')
2.
dsolve('Dy=x','y(0)=0','x')
1.MATLAB基本操作：矩阵Fra bibliotek乘、矩阵LU分解及矩阵分解为另一个矩阵和其转置相乘、解矛盾方程、符号求解微分方程。
（1）矩阵相乘：A=[3 4;5 6] B=[0 8;1 6]，求A*B。
（2）矩阵分解：A=[2 3;4 5]，求[L U]=lu(A)、B=Chol(A)。
（3）解方程：A=[2 8;3 4] b=[6;2]，求X=A\b。
ans =
-x-1+exp(x)*C1
>> dsolve('Dy=x','x')
ans =
1/2*x^2+C1
实验报告七常微分方程初值问题的数值解法
浙江大学城市学院实验报告课程名称数值计算方法实验项目名称常微分方程初值问题的数值解法 实验成绩指导老师签名日期2015/12/16 一.实验目的和要求1． 用Matlab 软件掌握求微分方程数值解的欧拉方法和龙格－库塔方法； 2． 通过实例学习用微分方程模型解决简化的实际问题;二.实验内容和原理编程题2-1要求写出Matlab 源程序m 文件,并有适当的注释语句；分析应用题2-2,2-3,2-4,2-5要求将问题的分析过程、Matlab 源程序和运行结果和结果的解释、算法的分析写在实验报告上; 2-1 编程编写用向前欧拉公式和改进欧拉公式求微分方程数值解的Matlab 程序,问题如下：在区间[],a b 内(1)N +个等距点处,逼近下列初值问题的解,并对程序的每一句添上注释语句; Euler 法y=eulera,b,n,y0,f,f1,b1改进Euler 法y=eulerproa,b,n,y0,f,f1,b1 2-2 分析应用题假设等分区间数100n =,用欧拉法和改进欧拉法在区间[0,10]t ∈内求解初值问题()()20(0)10y t y t y '=-⎧⎨=⎩并作出解的曲线图形,同时将方程的解析解也画在同一张图上,并作比较,分析这两种方法的精度; 2-3 分析应用题用以下三种不同的方法求下述微分方程的数值解,取10h = 画出解的图形,与精确值比较并进行分析; 1欧拉法； 2改进欧拉法； 3龙格－库塔方法；2-4 分析应用题考虑一个涉及到社会上与众不同的人的繁衍问题模型;假设在时刻t 单位为年,社会上有人口()x t 人,又假设所有与众不同的人与别的与众不同的人结婚后所生后代也是与众不同的人;而固定比例为r 的所有其他的后代也是与众不同的人;如果对所有人来说出生率假定为常数b ,又如果普通的人和与众不同的人的婚配是任意的,则此问题可以用微分方程表示为：其中变量()()()i p t x t x t =表示在时刻t 社会上与众不同的人的比例,()i x t 表示在时刻t 人口中与众不同的人的数量;1假定(0)0.01,0.02p b ==和0.1r =,当步长为1h =年时,求从0t =到50t =解()p t 的近似值,并作出近似解的曲线图形;2精确求出微分方程的解()p t ,并将你当50t =时在分题b 中得到的结果与此时的精确值进行比较; MATLAB 相关函数求微分方程的解析解及其数值的代入dsolve‘egn1’,‘egn2’,‘x ’ subsexpr,{x,y,…},{x1,y1,…}其中‘egn i ’表示第i 个方程,‘x ’表示微分方程中的自变量,默认时自变量为t ; subs 命令中的expr 、x 、y 为符合型表达式,x 、y 分别用数值x1、x2代入; >>symsxyz>>subs'x+y+z',{x,y,z},{1,2,3} ans= 6>>symsx>>subs'x^2',x,2 ans= 4>>s=dsolve‘12Dy y ∧=+’,‘(0)1y =’,‘x ’ ans= >>symsx >>subss,x,2 ans=右端函数(,)f x y 的自动生成f=inline ‘expr ’,’var1’,‘var2’,……其中’expr ’表示函数的表达式,’var1’,‘var2’表示函数表达式中的变量,运行该函数,生成一个新的函数表达式为fvar1,var2,……; >>f=inline'x+3y','x','y' f=Inlinefunction: fx,y=x+3y >>f2,3 ans= 114,5阶龙格－库塔方法求解微分方程数值解t,x=ode45f,ts,x0,options其中f 是由待解方程写成的m 文件名；x0为函数的初值；t,x 分别为输出的自变量和函数值列向量,t的步长是程序根据误差限自动选定的;若ts=t0,t1,t2,…,tf,则输出在自变量指定值,等步长时用ts=t0:k:tf,输出在等分点；options 用于设定误差限可以缺省,缺省时设定为相对误差310-,绝对误差610-,程序为：options=odeset ‘reltol ’,rt,’abstol ’,at,这里rt,at 分别为设定的相对误差和绝对误差;常用选项见下表;选项名 功能 可选值 省缺值 AbsTol 设定绝对误差正数 RelTol 设定相对误差 正数InitialStep 设定初始步长 正数 自动 MaxStep设定步长上界正数MaxOrder 设定ode15s 的最高阶数 1,2,3,4,5 5 Stats 显示计算成本统计 on,off off BDF 设定ode15s 是否用反向差分on,offoff例：在命令窗口执行>>odefun =inline ‘2*y t y -’,‘t ’,‘y ’；>>[],45(,[0,4],1)t y ode odefun =；ans=>>t y ‘o-’,%解函数图形表示>>45(,[0,4],1)ode odefun %不用输出变量,则直接输出图形 >>[],45(,0:4,1)t y ode odefun =；[],t yans=三.操作方法与实验步骤包括实验数据记录和处理2-1编程编写用向前欧拉公式和改进欧拉公式求微分方程数值解的Matlab 程序,问题如下：在区间[],a b 内(1)N +个等距点处,逼近下列初值问题的解,并对程序的每一句添上注释语句; Euler 法y=eulera,b,n,y0,f,f1,b1改进Euler 法y=eulerproa,b,n,y0,f,f1,b1Euler 法y=eulera,b,n,y0,f,f1,b1 y=zeros1,n+1; y1=y0; h=b-a/n; x=a:h:b; fori=1:n; yi+1=yi+hfxi,yi; end plotx,y holdon%求微分方程的精确解 x1=linspacea,b,100; '精确解为' s=dsolvef1,b1,'x' symsxy1=zeros1,100; for i=1:100y1i=subss,x,x1i; endplotx1,y1,'r'title'红色代表精确解'改进Euler 法y=eulerproa,b,n,y0,f,f1,b1 %求微分方程的数值解 y=zeros1,n+1; y1=y0; h=b-a/n; x=a:h:b; fori=1:n; T1=fxi,yi; T2=fxi+1,yi+hT1; yi+1=yi+h/2T1+T2; end plotx,y holdon%求微分方程的精确解 x1=linspacea,b,100; '精确解为' s=dsolvef1,b1,'x' symsxy1=zeros1,100; fori=1:100 y1i=subss,x,x1i; endplotx1,y1,'r'title'红色代表精确解' 2-2分析应用题假设等分区间数100n =,用欧拉法和改进欧拉法在区间[0,10]t ∈内求解初值问题()()20(0)10y t y t y '=-⎧⎨=⎩并作出解的曲线图形,同时将方程的解析解也画在同一张图上,并作比较,分析这两种方法的精度;1向前欧拉法>>euler0,10,100,10,inline'y-20','x','y','Dy=y-20','y0=10' ans= 精确解为 s= 20-10expx ans= +005Columns1through8（2）改进欧拉法>>eulerpro0,10,100,10,inline'y-20','x','y','Dy=y-20','y0=10' ans= 精确解为 s= 20-10expx ans= +005Columns1through8改进欧拉法的精度比向前欧拉法更高; 2-3分析应用题用以下三种不同的方法求下述微分方程的数值解,取10h = 画出解的图形,与精确值比较并进行分析; 1欧拉法； 2改进欧拉法；2-4分析应用题考虑一个涉及到社会上与众不同的人的繁衍问题模型;假设在时刻t 单位为年,社会上有人口()x t 人,又假设所有与众不同的人与别的与众不同的人结婚后所生后代也是与众不同的人;而固定比例为r 的所有其他的后代也是与众不同的人;如果对所有人来说出生率假定为常数b ,又如果普通的人和与众不同的人的婚配是任意的,则此问题可以用微分方程表示为：其中变量()()()i p t x t x t =表示在时刻t 社会上与众不同的人的比例,()i x t 表示在时刻t 人口中与众不同的人的数量;1假定(0)0.01,0.02p b ==和0.1r =,当步长为1h =年时,求从0t =到50t =解()p t 的近似值,并作出近似解的曲线图形;2精确求出微分方程的解()p t ,并将你当50t =时在分题b 中得到的结果与此时的精确值进行比较;1>>euler0,50,50,,inline'','t','p','Dp=','p0= 1' ans= 精确解为 s=1-99/100expx/500 ans=Columns1through82>>dsolve'Dp=','p0=','t' ans=1-99/100expt/500 >>1-99/100exp ans=与欧拉法求得的精确值差0,0001四.实验结果与分析。
微分方程数值解实验报告
微分方程数值解实验报告微分方程数值解实验报告一、引言微分方程是数学中一类重要的方程，广泛应用于各个科学领域。
在实际问题中，往往难以得到微分方程的解析解，因此需要借助数值方法来求解。
本实验旨在通过数值解法，探索微分方程的数值解及其应用。
二、数值解法介绍常用的微分方程数值解法有欧拉法、改进欧拉法、四阶龙格-库塔法等。
在本实验中，我们将采用改进欧拉法进行数值解的求取。
改进欧拉法是一种一阶的显式迭代法，其基本思想是将微分方程的导数用差商来近似表示，并通过迭代逼近真实解。
具体迭代公式如下：\[y_{n+1} = y_n + h \cdot f(x_n + \frac{h}{2}, y_n + \frac{h}{2} \cdot f(x_n, y_n))\]其中，\(y_n\)表示第n步的近似解，\(h\)表示步长，\(f(x_n, y_n)\)表示微分方程的导数。
三、实验步骤1. 确定微分方程及初始条件在本实验中，我们选择经典的一阶常微分方程：\[y' = -2xy\]并给定初始条件\(y(0) = 1\)。
2. 设置步长和迭代次数为了得到较为准确的数值解，我们需要合理选择步长和迭代次数。
在本实验中，我们将步长设置为0.1，迭代次数为10。
3. 迭代计算数值解根据改进欧拉法的迭代公式，我们可以通过编写计算程序来求解微分方程的数值解。
具体计算过程如下：- 初始化：设定初始条件\(y_0 = 1\)，并给定步长\(h = 0.1\)。
- 迭代计算：使用改进欧拉法的迭代公式，通过循环计算得到\(y_1, y_2, ...,y_{10}\)。
- 输出结果：将计算得到的数值解输出，并进行可视化展示。
四、实验结果与分析通过以上步骤，我们得到了微分方程的数值解\(y_1, y_2, ..., y_{10}\)。
将这些数值解进行可视化展示，可以更直观地观察到解的变化趋势。
根据实验结果，我们可以发现随着迭代次数的增加，数值解逐渐逼近了真实解。
第5次实验报告(常微分方程初值问题的数值解法)
班级： 学号： 姓名： 成绩：实验5 常微分方程初值问题的数值解法实验1实验目的1）熟悉欧拉法、改进欧拉法和龙格-库塔法的原理。
2）根据以上方法，编程求解常微分方程初值问题。
2 实验内容（1）编写程序，用以上各种方法求解教材P232例7-1、习题6、11的初值问题。
(2) 使用系统自带的函数dsolve 和ode45求例7-1的符号解析解和数值解。
3实验原理求解微分方程初值问题00(,)()y f x y y x y '=⎧⎨=⎩ （1） 欧拉法（显式）：10(,)n n n n n y y hf x y x x nh +=+⎧⎨=+⎩（2） 改进欧拉法(0)1(0)111(,)(,)(,)2n n n n n n n n n n y y hf x y h y y f x y f x y ++++⎧=+⎪⎨⎡⎤=++⎪⎣⎦⎩ （3） 经典龙格-库塔法（四阶）11234121324300(22)6(,)(,)22(,)22(,)()i i i i i i i i i i h y y K K K K K f x y h h K f x y K h h K f x y K K f x h y hK y y x +⎧=++++⎪⎪=⎪⎪=++⎪⎨⎪=++⎪⎪⎪=++⎪=⎪⎩4实验步骤1)建立函数文件，根据各公式编写程序；2)上机调试程序，运行程序进行计算，记录计算结果；3)分析各公式计算结果，比较各公式的优缺点。
5 程序设计欧拉法改进欧拉法function Euler1(x0,y0,h,n)%(x0,y0)：方程的初值%h：步长%n：计算的步数for i=1:nx=x0+h;yp=y0+h*f(x0,y0);yc=y0+h*f(x,yp);y=(yp+yc)/2;x %在屏幕显示每一步的x值y %在屏幕显示每一步计算的方程的数值解 x0=x;y0=y;end经典龙格-库塔法1）函数function f=f(x,y) f=y-2*x/y;end6总结注：若要更改matlab计算的数值类型，可以通过在matlab中设置实现：File -> Preferences ->Array Editor窗口中，Format 下方将Default array format设置为：long解微分方程的MATLAB命令MATLAB中主要用dsolve求微分方程的符号解析解，ode45求数值解。
数值分析实验报告之常微分方程数值解
数学与计算科学学院实验报告
实验项目名称常微分方程数值解
所属课程名称数值方法B
实验类型验证
实验日期2013.11.11
班级
学号
姓名
成绩
图1 h=0.1时三个方法走势图
图2 h=0.05时三个方法走势图
图4 h=0.05时三个方法走势图
附录1：源程序
附录2：实验报告填写说明
1．实验项目名称：要求与实验教学大纲一致。
2．实验目的：目的要明确，要抓住重点，符合实验教学大纲要求。
3．实验原理：简要说明本实验项目所涉及的理论知识。
4．实验环境：实验用的软、硬件环境。
5．实验方案（思路、步骤和方法等）：这是实验报告极其重要的内容。
概括整个实验过程。
对于验证性实验，要写明依据何种原理、操作方法进行实验，要写明需要经过哪几个步骤来实现其操作。
对于设计性和综合性实验，在上述内容基础上还应该画出流程图、设
计思路和设计方法，再配以相应的文字说明。
对于创新性实验，还应注明其创新点、特色。
6．实验过程（实验中涉及的记录、数据、分析）：写明具体实验方案的具体实施步骤，包括实验过程中的记录、数据和相应的分析。
7．实验结论（结果）：根据实验过程中得到的结果，做出结论。
8．实验小结：本次实验心得体会、思考和建议。
9．指导教师评语及成绩：指导教师依据学生的实际报告内容，给出本次实验报告的评价。
数值分析实验报告之常微分方程数值解
数学与计算科学学院告验报实常微分方程数值解实验项目名称数值方法B 所属课程名称验证实验类型2013.11.11 期验实日级班学号姓名成绩专业资料..资料．..资料．..资料．..资料．.1.395322 1.306570 0.9 1.348678 1.3072281.4581271.01.3685411.4138111.367880时三个方法走势图1 h=0.1图）i的范围进行扩大，即for(i=0;i<20;i++)h=0.05（此时将源程序中时三个方法与精确值的真值表2 h=0.05表精确值经典四阶库Euler法预估校正法步长1.011721 1.001229 0.05 1.001250 1.002500 1.024908 0.10 1.004837 1.004877 1.007375 1.039504 0.15 1.010708 1.014506 1.010764 1.055455 0.20 1.018731 1.023781 1.018802 1.072710 1.028801 0.25 1.028885 1.035092 1.0912171.0408181.0483370.301.040915.资料．..资料．.h=0.05时三个方法走势图图2dy?x?yxe??? 2、dx??1(0)?y?1由题可知精确解为：x?2e?(x2)。
x=0时，y(x)=0，当2h=0.1时三个方法与精确值的真值表3 h=0.1表精确值经典四阶库Euler 法预估校正法步长0.1 0.9295330.9094280.9000000.9096250.8725640.8355930.2 0.835927 0.8192490.7756550.3 0.754433 0.8268220.776081 0.727671 0.702726 0.7903480.4 0.727189.资料．.0.7614570.6881270.6886360.50.6617260.6572250.60.7387090.6567110.6293960.6314530.6319570.7208740.6040180.70.80.6111000.7069080.6115820.5841470.6959270.5945990.5950500.90.5685750.5814870.6871910.5562970.5810721.0时三个方法走势图h=0.1图3）的范围进行扩大，即for(i=0;i<20;i++)h=0.05（此时将源程序中i时三个方法与精确值的真值表h=0.05表4精确值预估校正法经典四阶库步长Euler法.资料．..资料．.图4 h=0.05时三个方法走势图【实验结论】（结果）1.预估校正法的精确度比经典的四阶库法的精确度高，库塔法最低；2.从表中数据可知三个算法所得数据与精确值相比，可得出以下结论（针对方案二）：（1）Euler法所得值偏离精确值最大，因此可知其精度相对来说最差；（2）经典的四阶库塔法所得值与精确值距离较近，因此可知对于Euler来说，该法更加有效；（3）预估校正法的数据时距离精确值最近的，其骤减幅度较小，因此对精度上的考虑而言，预估校正法应属于最佳解法；（4）由数据可知，上述两个方程的解的光滑性都比较差，从而导致四阶龙格库塔法的精度低于预估校正法的精度。
微分方程数值解法实验报告
微分方程数值解法实验报告实验题目：数值解微分方程的实验研究引言：微分方程是描述自然现象、科学问题和工程问题中变量之间的关系的重要数学工具。
然而，大部分微分方程很难找到解析解，因此需要使用数值方法来近似求解。
本实验旨在研究数值解微分方程的方法和工具，并通过实验验证其有效性和准确性。
实验步骤：1.了解微分方程的基本概念和求解方法，包括欧拉法、改进的欧拉法和龙格-库塔法。
2. 配置实验环境，准备实验所需的工具和软件，如Python编程语言和相关数值计算库。
3.选择一种微分方程进行研究和求解，可以是一阶、二阶或更高阶的微分方程。
4.使用欧拉法、改进的欧拉法和龙格-库塔法分别求解选定的微分方程，并比较其结果的准确性和稳定性。
5.计算数值解与解析解之间的误差，并进行误差分析和讨论。
6.对比不同数值解法的性能，包括计算时间和计算精度。
7.结果展示和总结，根据实验结果对数值解方法进行评价和选取。
实验结果：以一阶线性常微分方程为例，我们选择经典的“衰减振荡”问题进行实验研究。
该问题的微分方程形式为：dy/dt = -λy其中，λ为正实数。
我们首先使用Python编程语言实现了欧拉法、改进的欧拉法和龙格-库塔法。
进一步，我们选择了λ=0.5和初始条件y(0)=1，使用这三种数值解法求解该微分方程，并比较结果的准确性。
通过对比数值解和解析解可以发现，在短时间内，欧拉法、改进的欧拉法和龙格-库塔法的结果与解析解非常接近。
但随着时间的增加，欧拉法的结果开始偏离解析解，而改进的欧拉法和龙格-库塔法仍然能够提供准确的近似解。
这是因为欧拉法采用线性逼近的方式，误差随着步长的增加而累积，而改进的欧拉法和龙格-库塔法采用更高阶的逼近方式，可以减小误差。
为了更直观地比较不同方法的性能，我们还计算了它们的计算时间。
实验结果显示，欧拉法计算时间最短，而龙格-库塔法计算时间最长。
这表明在计算时间要求较高的情况下，可以选择欧拉法作为数值解方法。
微分方程数值解试验报告
中国矿业大学（北京）理学院微分方程数值解实验报告实验名称：常微分方程数值解法学号： XXX专业年级：信息与计算科学11级姓名： XXX实验时间： 2013年3月成绩：一、实验目的，内容 二、相关背景知识介绍 三、代码四、数值结果与分析 六、计算结果的分析一、实验目的，内容掌握一阶常微分方程的初值问题常用数值解法编写程序，利用Euler 法、改进Euler 法、Adams 外插内插法及Runge-Kutta 法解决相关问题。
二、相关背景知识介绍研究一阶常微分方程的初值问题的数值解1、Euler 法和改进Euler 法 （1）Euler 法用差商近似导数问题转化为（2）改进Euler 法]2/),(),([111++++=+n n n n n u t f u t f h y y n00(,)()dy f x y a x bdx y x y ⎧=⎪≤≤⎨⎪=⎩dxdy h y y n n ≈-+1,...)3,2,1,0()(),(001=⎩⎨⎧=+=+n x y y y x hf y y n n n n（3）线性多步法（Adams 外插内插法）j n k j j j n k j j f h u+=+=∑∑=00βα （4）Runge-Kutta 法用f 在一些点的值表示ψ（tn,un,h ）,使单步局部截断误差的阶和Taylor 展开法相等。
将初值问题写成积分形式：di i u i f t u h t u h t t ⎰++=+))(,()()(三、代码（Matlab ）题目：用Euler 法、改进Euler 法、Adams 外插内插法及Runge-Kutta 法解决下面一阶常微分方程初值问题。
代码:方程：function y=f(t,u)y=2*t*u/(1+t^2);endEuler 法:function y=euler(n,u0,x0)t=1;u=zeros(n+1,1);h=t/n;u(1)=u0;x(1)=x0;for i=1:nx(i+1)=x0+h*i;u(i+1)=u(i)+f(x(i),u(i))*h;endPlot(x,u,'r*')end改进Euler法:function y=imeuler(n,u0,x0)err=10^-6;t=1;u=zeros(n+1,1);x=zeros(n+1,1);h=t/n;u(1)=u0;x(1)=x0;for i=1:nx(i+1)=x0+h*i;u(i+1)=u(i);ut=30;while abs(u(i+1)-ut)>errut=u(i+1);u(i+1)=u(i)+(f(x(i),u(i))+f(x(i+1),u(i+1)))*0.5*h;endendplot(x,u)endAdams内插法：k=3function y=iadams(n,u0,x0)t=1;u=zeros(n+1,1);h=t/n;u(1)=u0;u(2)=u0;u(3)=u0;x(1)=x0;x(2)=x0+h;x(3)=x0+2*h;for i=3:nx(i+1)=x0+h*i;u(i+1)=u(i)+h*(9*f(x(i+1),u(i+1))+19*f(x(i),u(i))-5*f(x(i-1),u(i-1))+f(x(i-2),u(i-2)))/24;endplot(x,u)endRunge-Kutta法：（本文采用heun法）function y=heun(n,u0,x0)t=1;u=zeros(n+1,1);h=t/n;u(1)=u0;x(1)=x0;for i=1:nx(i+1)=x0+h*i;k1=f(x(i),u(i));k2=f(x(i)+0.3333*h,u(i)+0.33333*h*k1);k3=f(x(i)+0.666663*h,u(i)+0.666666*h*k2);u(i+1)=u(i)+0.25*h*(k1+3*k3);endplot(x,u)end四、数值结果与分析对于所给方程，各个方法均给出了良好的解，其实heun方法最好10步就出了的结果。
常微分方程的求解与定性分析实验报告
常微分方程的求解与定性分析实验报告 一、实验综述1、实验目的及要求● 归纳和学习求解常微分方程(组)的基本原理和方法；● 掌握解析、数值解法，并学会用图形观察解的形态和进行解的定性分析； ● 熟悉MATLAB 软件关于微分方程求解的各种命令；● 通过范例学习建立微分方程方面的数学模型以及求解全过程；● 通过该实验的学习，使学生掌握微分方程(组)求解方法（解析法、欧拉法、梯度法、改进欧拉法等），对常微分方程的数值解法有一个初步了解，同时学会使用MATLAB 软件求解微分方程的基本命令，学会建立微分方程方面的数学模型。
这对于学生深入理解微分、积分的数学概念，掌握数学的分析思维方法，熟悉处理大量的工程计算问题的方法是十分必要的。
2、实验仪器、设备或软件电脑 、matlab7.0二、实验过程（实验步骤、记录、数据、分析）实验内容：根据实验内容和步骤，完成以下实验，要求写出实验报告（实验目的→问题→数学模型→算法与编程→计算结果→分析、检验和结论）1．求微分方程的解析解，并画出它们的图形。
y '= y + 2 x, y (0) = 1, 0< x <1；m=dsolve('Dy=y+2*x','y(0)=1','x')ezplot(m,[0 1])m =3*exp(x) - 2*x – 21．求微分方程⎪⎩⎪⎨⎧====-+]100[0)0(;0)0(01.03t uu u u u &&& 的数值解，要求编写求解程序。
function dy=vdp1000(t,y)dy=zeros(2,1);dy(1)=y(2);dy(2)=-y(1)+0.1*y(1)^3;[T,Y]=ode15s('vdp1000',[0 10],[0 0]);plot(T,Y(:,1),'-')3．Rossler 微分方程组：当固定参数b =2，c =4时，试讨论随参数a 由小到大变化（如 a ∈(0，0.65))而方程解的变化情况，并且画出空间曲线图形，观察空间曲线是否形成混沌状？function r=rossler(t,x)global a;global b;global c;r=[-x(2)-x(3);x(1)+a*x(2);b+x(3)*(x(1)-c)];global a;global b;global c;b=2;c=4;t0=[0,200];for a=0:0.1:0.6[t,x]=ode45('rossler',t0,[0,0,0]);subplot(1,2,1);plot(t,x(:,1),'r',t,x(:,2),'g',t,x(:,3),'b');title('x(红色),y(绿色),z(蓝色)随t的变化情况');xlabel('t');subplot(1,2,2);plot3(x(:,1),x(:,2),x(:,3))title('相图');xlabel('x');ylabel('y');zlabel('z');pauseend结果显示：a=0:a=0.1:a=0.2:a=0.3:a=0.4:a=0.5:结果分析：从图像可以看出，当a=0时，微分方程的解（x,y,z）收敛与（0，0.5，0.5）；当a=0.1时，（x,y,z）仍收敛与（0，0.5，0.5），只是收敛速度减慢；当a=0.2时，（x,y,z）已发散，周期性变化；随着a的增大，(x,y,z)接近其极限环的速度加快，空间曲线成混沌状。
数学实验报告——利用MALTAB计算常微分方程数值解
实验二 常微分方程数值解一、火箭飞行器㈠问题描述小型火箭初始质量为1400kg ，其中包括1080kg 燃料，火箭竖直向上发射时燃料燃烧率为18kg/s ，由此产生32000N 的推力，火箭引擎在燃料用尽时关闭，设火箭上升时空气阻力正比于速度的平方，比例系数为0.4kg/m ，求引擎关闭瞬间火箭的高度，速度，加速度及火箭达到最高点时的高度和加速度，并画出高度、速度、加速度随时间变化的图形。
㈡方法与公式1、简要分析本题的求解需要用到常微分方程，而整个过程又被分为两个阶段：火箭加速上升阶段和燃料燃尽后减速的阶段。
由题目易知第一个阶段持续时间T 1=108018=60s 列出第一阶段的方程组：设M0为火箭本身质量，m 为燃料质量，g 为重力加速度 = 9.8m/s 2,燃料燃烧率为a ，空气阻力的比例系数为k ，F 为推进力。
M0 = 1400-1080 = 320kg; v̇=F−kv 2−(M0+m )g M0+mm =−a初值v̇=0，m =1080。
由以上各式可以求出t=T 1时火箭的速度。
再求解第二阶段：v̇=−kv 2−M0g M0m =0可以求出火箭速度降为0的时刻。
将整个过程中的时间向量以及速度向量联合起来，利用第三章所学插值与数值积分的方法可以求得任意时刻火箭的近似高度。
2、方法求解常微分方程时，我分别采用了自己编写的欧拉公式、改进欧拉公式、4级4阶龙格-库塔公式，以及MATLAB自带的龙格-库塔方法，求解数值积分时采用辛普森公式。
由于Matlab自带的Simpson公式是自适应的，因此需要使用自己在上一次实验时所编的Simpson公式。
㈢结果与分析1、各种公式的对比首先，我作出了各种不同公式计算得到的火箭速度随时间变化的图像，图如下：从图中可以看出，各种公式计算得到的结果基本一致，为确定其区别，将图像放大，放大约2000 倍后，得到下图：分析：从图中可以看出，自编欧拉公式距离MATLAB 自带龙格-库塔公式最远，精度最差；自编的改进欧拉公式和自编的龙格-库塔公式结果基本一致，两者中自编龙格-库塔公式距MATLAB 自带龙格-库塔公式的结果稍近。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
常微分方程的数值解法 专业班级：信息软件 姓名：吴中原 学号：120108010002 一、实验目的 1、熟悉各种初值问题的算法，编出算法程序； 2、明确各种算法的精度与所选步长有密切关系；通过计算更加了解各种 算法的优越性。
二、实验题目 1、根据初值问题数值算法，分别选择二个初值问题编程计算； 2、试分别取不同步长，考察某节点jx处数值解的误差变化情况； 3、试用不同算法求解某初值问题，结果有何异常； 4、分析各个算法的优缺点。
三、实验原理与理论基础
（一） 欧拉法算法设计 对常微分方程初始问题
(6-1) (6-2) 用数值方法求解时，我们总是认为(6-1)、(6-2)的解存在且唯一。 欧拉法是解初值问题的最简单的数值方法。从(6-2)式由于y (x0) = y0已给定，因而可以算出),()('000yxfxy。 设x1 = h充分小，则近似地有：
),()(')()(00001yxfxyhxyxy (6-3)
记 ,n,,ixyyii10 )(
从而我们可以取 ),(0001yxhfyy 作为)(1xy的近似值。利用1y及f (x1, y1)又可以算出)(2xy的近似值： ),(1112yxhfyy 一般地，在任意点hnxn11处)(xy的近似值由下式给出 ),(1nnnnyxhfyy (6-4)
这就是欧拉法的计算公式，h称为步长。
)( ),(dd00yxyyxf
x
y （二）四阶龙格-库塔法算法设计： 欧拉公式可以改写为：111 ,iiiiyykkhfxy，它每一步计算,fxy的值一次，截
断误差为2oh。
改进的欧拉公式可以改写为：11212112,,iiiiiiyykkkhfxykhfxhyk ，它每一步要计算,fxy的值两次，截断误差为3oh。
改进的欧拉方法之所以比欧拉方法具有更高的精度，是因为在每一步它都比欧拉方法多计算了一次,fxy的值。因此，要进一步提高精度，可以考虑在每一步增加计算,fxy的次数。 如果考虑在每一步计算,fxy的值四次，则可以推得如下公式：
1123412132431226, 1,221,22,iiiiiiiiiiyykkkkkhfxyhkhfxykhkhfxykkhfxhyk
此公式称为标准四阶龙格-库塔(Runge-Kutta)公式，它的截断误差为5oh。虽然用龙格-库塔方法每一步需要四次调用f，计算量较改进的欧拉方法大一倍，这里由于龙格-库塔方法的步长增大了一倍，因而两种方法总的计算量相同，但龙格-库塔方法精确度更高。所以龙格-库塔公式兼顾了精度和计算工作量的较为理想的公式，在实际计算中最为常用。 四、实验内容
（一）问题重述： 科学计算中经常遇到微分方程（组）初值问题，需要利用Euler法，改进Euler 法，Rung-Kutta方法求其数值解，诸如以下问题： （1）0 y xy yxy 04 0 < x1 分别取h=0.1,0.2,0.4时数值解。
初值问题的精确解254xey
。
（2）用r=3的Adams显式和预 - 校式求解
2201yxy
y 01x
取步长h=0.1，用四阶标准R-K方法求值。 （3）用改进Euler法或四阶标准R-K方法求解
10330012101213y yy 2y yy y y y
10x 取步长h0.01,计算15.0,10.0,05.0y y y数值解，参考结果
8613125.015.0,1493359.015.0,9880787.015.0321y y y
（4）利用四阶标准R- K方法求二阶方程初值问题的数值解
（I）
10,00023yyyyy
0.02h x10
(II) 00,100)1(1.02y yyyyy 0.1h x,10 (III) 00,101yyeyyx 0.1h x,20 (IV) 00,100sinyyyy 2,400.h x (二)实验代码：
1、欧拉法程序
function y=Euler(a,b,M,y0) %a=1,b=2,M=10,f=t*y^(1/3),y0=1; h=(b-a)/M; t=zeros(1,M+1); t=a:h:b; y=zeros(1,M+1); yy=zeros(1,M+1); y(1)=y0; for k=1:M y(k+1)=y(k)+h*t(k)*y(k)^(1/3); end yb=y(M+1); yy=((t.^2+2)./3).^1.5; det=yy-y; plot(t,y,'r-',t,yy,'b:',t,det); 2、改进欧拉法程序 function H=heeuler(a,b,M,ya,f) %a=0,b=1,M=10,f=t*t+t-y,y0=0; h=(b-a)/M; t=zeros(1,M+1); y=zeros(1,M+1); p=0;q=0; t=a:h:b; y(1)=ya; for k=1:M p=feval(f,t(k),y(k)); q=feval(f,t(k+1),y(k)+h*p); y(k+1)=y(k)+0.5*h*(p+q); end yy=t.*t-t+1-exp(-t); det=yy-y; plot(t,y,'r-',t,yy,'b:',t,det); H=[t',y',yy',det']
function f=ff(t,y); f=t.^2+t-y; 3、四阶龙格-库塔法程序 function H=r_k4(a,b,M,ya,f) %a=0,b=1,M=10,f=t*t+t-y,y0=0; h=(b-a)/M; t=zeros(1,M+1); t=a:h:b; y=zeros(1,M+1); K1=0;K2=0;K3=0;K4=0; y(1)=ya; for k=1:M K1=feval(f,t(k),y(k)); K2=feval(f,t(k)+0.5*h,y(k)+0.5*h*K1); K3=feval(f,t(k)+0.5*h,y(k)+0.5*h*K2); K4=feval(f,t(k)+h,y(k)+h*K3); y(k+1)=y(k)+1/6*(K1+2*K2+2*K3+K4); end yy=t.*t-t+1-exp(-t); det=yy-y; plot(t,y,t,yy,t,det); H=[t',y',yy',det'] function f=ff(t,y); f=t.^2+t-y;
五、实验结果
1）0 y xy yxy 04 0 < x1
分别取h=0.1,0.2,0.4时数值解。
初值问题的精确解254xey
。
Euler('han',0,0,2,0.1) ans = 0 0 0.1000 0 0.2000 0.0400 0.3000 0.1192 0.4000 0.2356 0.5000 0.3862 0.6000 0.5669 0.7000 0.7729 0.8000 0.9988 0.9000 1.2389 1.0000 1.4874 1.1000 1.7386 1.2000 1.9874 1.3000 2.2289 1.4000 2.4591 1.5000 2.6749 1.6000 2.8736 1.7000 3.0539 1.8000 3.2147 1.9000 3.3561 2.0000 3.4784 Euler('han',0,0,2,0.2) ans =
0 0 0.2000 0 0.4000 0.1600 0.6000 0.4672 0.8000 0.8911 1.0000 1.3886