常微分方程初值问题的仿真解法

常微分方程初值问题的解法及应用常微分方程是数学中非常重要的一部分，它涉及了许多领域的模型建立和问题求解。

本文将介绍常微分方程初值问题的解法及其应用。

一、常微分方程初值问题的定义常微分方程初值问题是指给定一个常微分方程，以及它在某一点上的初始条件，求解该方程的解曲线。

通常，一个常微分方程初值问题可以表示为:y'(x) = f(x,y), y(x0) = y0,其中，y(x)是未知函数，f(x,y)是已知函数，y(x0) = y0是初始条件。

二、常微分方程初值问题的解法常微分方程初值问题的解法有多种，下面我们将介绍几种常用的方法。

1.欧拉法欧拉法是最简单的一种求解常微分方程初值问题的方法。

该方法基于初始条件，通过不断迭代计算得到近似解曲线。

具体步骤如下：步骤1：设定步长h，确定求解区间[x0, xn]，计算步数n。

步骤2：初始化，即确定初始点(x0, y0)。

步骤3：根据方程dy/dx = f(x,y)和初始点(x0, y0)，计算斜率k = f(x0, y0)。

步骤4：根据已知的斜率和步长h，计算下一个点的坐标(xi+1,yi+1)。

步骤5：重复步骤3和步骤4，直到达到步数n。

步骤6：得到近似解曲线。

2.改进的欧拉法（改进欧拉法）改进的欧拉法是对欧拉法的改进，其求解精度比欧拉法更高。

具体步骤如下：步骤1：设定步长h，确定求解区间[x0, xn]，计算步数n。

步骤2：初始化，即确定初始点(x0, y0)。

步骤3：根据方程dy/dx = f(x,y)和初始点(x0, y0)，计算斜率k1 =f(x0, y0)。

步骤4：根据已知的斜率k1和步长h/2，计算中间点的坐标(x0+h/2, y0+k1*h/2)。

步骤5：根据方程dy/dx = f(x,y)和中间点的坐标(x0+h/2, y0+k1*h/2)，计算斜率k2= f(x0+h/2, y0+k1*h/2)。

步骤6：根据已知的斜率k2和步长h，计算下一个点的坐标(xi+1,yi+1)。

常微分方程初值问题的解法随着科技的不断进步和人类社会的不断发展，工程技术和科学技术的发展已经成为推动社会进步的重要力量，而数学则是工程技术和科学技术的基础和支撑，常微分方程作为数学分支的重要组成部分，对于理论研究和实际应用都有着深远的影响。

在实际工程中，解决常微分方程初值问题是数学理论在抽象式运算与工程实践之间的重要桥梁。

本文将介绍常微分方程初值问题的概念、求解方法以及实际应用。

一、常微分方程初值问题的概念常微分方程是指未知函数一阶或高阶微商与自变量和常数的关系式，常微分方程初值问题是指在初值u(x0)=u0已知的情况下，确定函数u(x)的解的问题。

在初值问题中，自变量是独立变量，取值范围可以是任意实数，因变量是函数值，是依赖自变量而实现的数值，常数是影响函数变化的一些固定参数。

常微分方程模型经常出现在工程技术模型中，一些实际应用场景可以通过建立数学模型来进行求解。

二、常微分方程初值问题的解法常微分方程初值问题的解法大致可以分为两种，一种是解析解法，即直接利用微积分学知识对方程进行求解；另一种是数值解法，即采用数值方法对方程进行数值计算求解。

下面将分别介绍这两种方法的解法原理。

1. 解析解法解析解法是指通过数学工具对函数解析表达式进行研究，以求出常微分方程的解。

该方法的先决条件是对方程具有严格的内部结构和特殊的形式，只有在特殊情况下才能找到一些特解。

这种方法的难点在于方程方程形式和初始条件可能存在巨大的数学难度，解析解的求解需要求解一些解析式的积分、微分和级数。

往往只有在一些特殊情况下，解析解法才能一般性的解决问题，因此该方法的适用场景相对较少。

2. 数值解法数值解法是指通过数值计算的方法，通过有限个代数运算和计算机模拟的方法得出方程的解。

数值解法的优点是具有广泛的适用性，可以有效地求解各种类型的常微分方程初值问题，使得无法通过解析方法求解的问题也可以得到解答。

数值解法可分为无条件稳定和条件稳定两种情况，前者是指方法不会出现不稳定结果的情况，而后者则保证了方法收敛性的同时，存在一定的条件限制。

当 x x1 时，得 y1 y0 f(x0 , y0 )(x1 x0 )
这么就取得了P1点旳坐标： (x1, y1) 。将y1作为y(x1)旳 近似值(想象(x1, y1) 在积分曲线y=y(x)上)
过点P1(x1,y1)，作积分曲线y=y(x)旳切线交直线x=x2于
P2点。注意切线 P1P2 旳斜率(近似)为 y(x1 ) f(x1 , y1 )
第9章 常微分方程初值问题数值解法
§9.1 引言
➢ 包括自变量、未知函数及未知函数旳导数旳方程称 为微分方程。
➢ 自变量个数只有一种旳微分方程称为常微分方 程。
微分方程中出现旳未知函数最高阶导数旳阶数 称为微分方程旳阶数。
假如未知函数y及其各阶导数
y, y, … , y(n)
都是一次旳,则称其为线性旳,不然称为非线性旳。
xi1 xi1 f[xi , y(xi )]
代入上式,并用yi近似替代式中y(xi)即可得到 两步欧拉公式
yi1 yi1 2hf(xi , yi ) ( 9.7 )
【注】欧拉措施和梯形措施，都是单步法,其特点是 在计算yi+1时只用到前一步旳信息yi; 而两步欧拉公式 (9.7)中除了yi外,还用到更前一步旳 信息yi-1,即调用了前两步旳信息。
当 x xi1 时,得
yi1 yi f(xi , yi )(xi1 xi )
这么,从x0逐一算出 x1 , x2 , … xn
相应旳数值解
y1 , y 2 , … yn
就取得了一系列旳点: P1, P1,…,Pn。 从图形上看,就取得了一条近似于曲线y=y(x)
旳折线 P1P2P3 … Pn 。
xi x0 ih, i 1,2, … , n
数值解法需要把连续性旳问题加以离散化，从 而求出离散节点旳数值解。

浙江大学城市学院实验报告课程名称数值计算方法实验项目名称常微分方程初值问题的数值解法 实验成绩指导老师签名日期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四.实验结果与分析。