数学实验第二次作业任务常微分方程数值求解
常微分方程数值解实验
有大量的常微分方程虽然从理论上讲，其解是存在的，但我们却无 法求出其解析解，此时，我们需要寻求方程的数值解，在求常微分方程 数值解方面，MATLAB具有丰富的函数，我们将其统称为solver，其一般 格式为：
Image
Image
如果微 分方程 由一个 或多个 高阶常微分方程给出，要得到该方程的数值解,可以将方程转换成一阶 常微分方程组。假设高阶常微分方程的一般形式为y( n) = f ( t, y, yʹ, ⋯,y( n - 1) )，而且函数y(t)的各阶导数初值为y(0)，yʹ(0) ,…, y( n - 1) (0)可以选 择一组变量令: x1= y, x2 = yʹ,…, xn = y( n - 1) ,我们就可以把原高阶常微 分方程转换成下面的一阶常微分方程组形式： 而且初值x1(0)=y(0),x2(0)=yʹ(0),…，xn(0)=(0)。 转换以后就可以求原 高阶常微分方程的数值解了。 例2 求微分方程，,的数值解。 对方程进行变换,选择变量 (1) 建立自定义函数 用edit命令建立自定义函数名为f.m,内容为: function y =f(t,x) y=[x(2);x(3);-t^2*x(2)*x(1)^2-t*x(1)*x(3)+exp(t*x(1))]; (2)调用对微分方程数值解ode45函数求解 用edit命令建立一个命令文件f2. m,内容为: >>x0=[2;0;0]; >>[t,y] =ode45(’f’,[0,10],x0);plot(t,y); >>figure; >>plot3(y(:,1),y(:,2), y(:,3))得
常微分方程的数值解法实验报告
常微分方程的数值解法专业班级：信息软件 姓名：吴中原 学号：120108010002 一、实验目的1、熟悉各种初值问题的算法，编出算法程序；2、明确各种算法的精度与所选步长有密切关系；通过计算更加了解各种 算法的优越性。
二、实验题目1、根据初值问题数值算法，分别选择二个初值问题编程计算；2、试分别取不同步长，考察某节点j x处数值解的误差变化情况； 3、试用不同算法求解某初值问题，结果有何异常； 4、分析各个算法的优缺点。
三、实验原理与理论基础（一） 欧拉法算法设计对常微分方程初始问题(6-1)(6-2)用数值方法求解时，我们总是认为(6-1)、(6-2)的解存在且唯一。
欧拉法是解初值问题的最简单的数值方法。
从(6-2)式由于y (x 0) = y 0已给定，因而可以算出),()('000y x f x y =。
设x 1 = h 充分小，则近似地有：),()(')()(00001y x f x y hx y x y =≈-(6-3)记 ,n ,,i x y y i i 10 )(== 从而我们可以取),(0001y x hf y y ==作为)(1x y 的近似值。
利用1y 及f (x 1, y 1)又可以算出)(2x y 的近似值：),(1112y x hf y y +=一般地，在任意点()h n x n 11+=+处)(x y 的近似值由下式给出),(1n n n n y x hf y y +=+(6-4)这就是欧拉法的计算公式，h 称为步长。
⎪⎩⎪⎨⎧==)( ),(d d 00y x y y x f x y（二）四阶龙格-库塔法算法设计：欧拉公式可以改写为：()111,i i i i y y k k hf x y +=+⎧⎪⎨=⎪⎩，它每一步计算(),f x y 的值一次，截断误差为()2o h 。
改进的欧拉公式可以改写为：()()()11212112,,i i i i i i y y k k k hf x y k hf x h y k +⎧=++⎪⎪=⎨⎪=++⎪⎩，它每一步要计算(),f x y 的值两次，截断误差为()3o h 。
实验报告——常微分方程的数值解法
实验报告实验项目名称常微分方程的数值解法实验室数学实验室所属课程名称微分方程数值解实验类型上机实验实验日期2013年3月11日班级10信息与计算科学学号2010119421姓名叶达伟成绩实验概述：【实验目的及要求】运用不同的数值解法来求解具体问题，并通过具体实例来分析比较各种常微分方程的数值解法的精度，为以后求解一般的常微分方程起到借鉴意义。
【实验原理】各种常微分方程的数值解法的原理，包括Euler法，改进Euler法，梯形法，Runge-Kutta方法，线性多步方法等。
【实验环境】（使用的软硬件）Matlab软件实验内容：【实验方案设计】我们分别运用Euler法，改进Euler法，RK方法和Adams隐式方法对同一问题进行求解，将数值解和解析解画在同一图像中，比较数值解的精度大小，得出结论。
【实验过程】（实验步骤、记录、数据、分析）我们首先来回顾一下原题：对于给定初值问题：1. 求出其解析解并用Matlab画出其图形；2. 采用Euler法取步长为0.5和0.25数值求解（2.16）,并将结果画在同一幅图中，比较两者精度;3. 采用改进Euler法求解（2.16），步长取为0.5；4. 采用四级Runge-Kutta法求解（2.16），步长取为0.5；5. 采用Adams四阶隐格式计算（2.16），初值可由四级Runge-Kutta格式确定。
下面，我们分五个步骤来完成这个问题：步骤一，求出（2.16）式的解析解并用Matlab 画出其图形； ，用Matlab 做出函数在上的图像，见下图：00.51 1.52 2.53 3.54 4.550.511.522.533.5x 1015y=exp(1/3 t 3-1.2t)exact solution图一 初值问题的解析解的图像步骤二，采用Euler 法取步长为0.5和0.25数值求解(2.16),并将结果画在同一幅图中，比较两者精度;我们采用Euler 法取步长为0.5和0.25数值求解，并且将数值解与解析解在一个图中呈现，见下图：00.51 1.52 2.53 3.54 4.550.511.522.533.5x 1015Numerical solution of Euler and exact solutionexact solution h=0.5h=0.25图二 Euler 方法的计算结果与解析解的比较从图像中不难看出，采用Euler 法取步长为0.5和0.25数值求解的误差不尽相同，也就是两种方法的计算精度不同，不妨将两者的绝对误差作图，可以使两种方法的精度更加直观化，见下图：00.51 1.52 2.53 3.54 4.550.511.522.533.5x 1015Absolute error of numerical solution and exact solutionh=0.5h=0.25图三 不同步长的Euler 法的计算结果与解析解的绝对误差的比较 从图像中我们不难看出，步长为0.25的Euler 法比步长为0.5的Euler 法的精度更高。
常微分方程数值解实验报告
常微分方程数值解实验报告实验报告：常微分方程数值解1.引言常微分方程（Ordinary Differential Equations, ODEs）是数学领域中一个重要的研究对象，涉及到许多自然科学和工程技术领域的问题。
解常微分方程的数值方法是一种求解差分方程的方法，通过计算机找到方程的近似解，对于模拟和预测连续过程非常有用。
本实验旨在通过数值解法，验证和应用常微分方程的解，并比较不同数值方法的精度和效率。
2.实验目的2.1理解常微分方程的基本概念和数值解法；2.2掌握将常微分方程转化为数值求解问题的基本方法；2.3运用数值解法求解常微分方程；2.4比较不同数值解法的精度和效率。
3.实验内容3.1 欧拉方法（Euler Method）给定一个一阶常微分方程dy/dx=f(x,y)，通过将其离散为差分形式，欧拉方法可以通过以下递推公式来求解：y_{n+1}=y_n+h*f(x_n,y_n)其中，h为步长，x_n和y_n为当前的x和y值。
3.2 改进的欧拉方法（Improved Euler Method）改进的欧拉方法使用欧拉方法的斜率的平均值来估计每一步中的斜率。
具体公式如下：k1=f(x_n,y_n)k2=f(x_n+h,y_n+h*k1)y_{n+1}=y_n+h*((k1+k2)/2)3.3 二阶龙格-库塔法（Second-order Runge-Kutta Method）二阶龙格-库塔法通过计算每个步骤中的两个斜率来估计每个步长中的斜率。
具体公式如下：k1=f(x_n,y_n)k2=f(x_n+h/2,y_n+(h/2)*k1)y_{n+1}=y_n+h*k24.实验步骤4.1选取常微分方程，并将其转化为数值求解问题的形式；4.2根据给定的初始条件和步长，使用欧拉方法、改进的欧拉方法和二阶龙格-库塔法求解该方程；4.3比较三种方法的数值解与理论解的差异，并分析其精度和效率；4.4尝试不同的步长，观察相应的数值解的变化。
常微分方程数值解实验
多步法，Gear’s反向
数值积分，精度中等
若ode45失效时，
可尝试使用
ode23s
刚性
一步法，2阶Rosebrock算法，
低精度。
当精度较低时，
计算时间比ode15s短
odefx为显式常微分方程 中的 ，t为求解区间，要获得问题在其他指定点 上的解，则令t=[t0,t1,t2,…]（要求 单调），y0初始条件。
MATLAB 中有几个专门用于求解常微分方程的函数,它们的设计思想基于Runge-Kutta方法,基本设计思想为：从改进的欧拉方法比欧拉方法精度高的缘由着手,如果在区间[ x1, xi+1]多取几个点的斜率值,然后求取平均值,则可以构造出精度更高的计算方法。 这些函数主要包括：ode45、ode23、ode15s、ode113、ode23s、ode23t、ode23tb. 其中最常用的是函数ode45,该函数采用变步长四阶五阶Runge-Kutta法求数值解,并采用自适应变步长的求解方法。ode23采用二阶三阶Runge-Kutta法求数值解,与ode45类似,只是精度低一些。ode15s用来求刚性方程组。
43
4月22日
588
666
28
46
4月23日
693
782
35
55
4月24日
774
863
39
64
4月25日
877
954
42
73
4月26日
988
1093
48
76
4月27日
1114
1255
56
78
4月28日
1199
1275
59
78
4月29日
常微分方程数值解法
而
y ( xn 1 ) y ( xn ) hy( xn ) 1 2 1 1 h y( xn ) h3 y ( xn ) y (4) ( xn ) O(h 5 ) 2! 3! 4!
要使y(xn+1)–yn+1= O(h5) ，比较h的同次幂系数，必 有
h : a 1 1 h : a b c d 1 2 h : a 2b 2d 1 h3 : a 3b 3d 1 h 4 : a 4b 4d 1
yn 1 yn hf ( xn , yn ) h yn 1 yn 2 [ f ( xn , yn ) f ( xn 1 , yn 1 )] (n 0,1, 2,)
局部截断误差为O(h3)，是二阶方法.
标准四阶Runge-Kutta公式
1 yn 1 yn ( K1 2 K 2 2 K 3 K 4 ) K hf ( x ,6y ) n n 1 h K1 ) K 2 hf ( xn , yn 2 2 K hf ( x h , y K 2 ) n n 3 2 2 K 4 hf ( xn h, yn K 3 )
用Euler预报—校正格式求解初值问题
dy y y 2 sin x 0 dx y (1.2)及y(1.4)的近似值， 小数点后至少保留5位. 解 设f(x,y)=-y-y2sinx , x0=1,y0=1, xn=x0+nh=1+0.2n, Euler预报—校正格式为
h yn 1 yn [ f ( xn , yn ) f ( xn1 , yn1 )] 2 0.2 yn 1 yn [8 3 yn 8 3 yn 1 ] 2
实验二 微分方程数值解
实验二微分方程数值解（2学时）
一．实验目的：
掌握欧拉法、四阶龙格库塔法的程序编制方法。
二．实验方法及预习内容：
对微分方程描述的控制系统，利用欧拉法、二阶龙格-库塔法、四阶龙格-库塔法分别编写M文件，进行数值计算和作图。
三．实验内容：
1．分别用欧拉法、二阶龙格-库塔法、四阶龙格-库塔法求下面系统的输出响应y(t)在0≤t≤1上，h=0.1时的数值解：
'2,(0)1
y y y
=-=
要求保留4位小数，并将三种方法的结果与真解2
()t
y t e-
=进行比较。
2. 四阶龙格-库塔法求下面系统的输出响应y(t)在2≤t≤3上，h=0.5时的数值解：
2
1(),23,(2)1
y t y t y
'=+-≤≤=
要求保留5位小数，并将结果与真解
1
1
y t
t
=+
-
进行比较。
四．实验总体要求：
1．每次实验应做好实验前的预习和准备；
2．实验后应及时提交仿真程序（M文件与Word文档）、实验结果和图示、实验分析与总结；
3．认真撰写实验报告。
五．本次实验要求：
熟悉微分方程数值解法（欧拉法、二阶与四阶龙格-库塔法）。
常微分方程的数值解法全文
第8章常微分方程的数值解法8.4单步法的收敛性与稳定性8.4.1相容性与收敛性上面所介绍的方法都是用离散化的方法，将微分方程初值问题化为差分方程初值问题求解的.这些转化是否合理？即当h →∞时，差分方程是否能无限逼近微分方程，差分方程的解n y 是否能无限逼近微分方程初值问题的准确解()n y x ，这就是相容性与收敛性问题.用单步法(8.3.14)求解初值问题(8.1.1)，即用差分方程初值问题100(,,)()n n n n y y h x y h y x y ϕ+=+⎧⎨=⎩(8.4.1)的解作为问题(8.1.1)的近似解，如果近似是合理的，则应有()()(,(),)0 (0)y x h y x x y x h h hϕ+--→→(8.4.2)其中()y x 为问题(8.1.1)的精确解.因为0()()lim ()(,)h y x h y x y x f x y h→+-'==故由(8.4.2)得lim (,,)(,)h x y h f x y ϕ→=如果增量函数(,(),)x y x h ϕ关于h 连续，则有(,,0)(,)x y f x y ϕ=(8.4.3)定义8.3如果单步法的增量函数(,,)x y h ϕ满足条件(8.4.3)，则称单步法(8.3.14)与初值问题(8.1.1)相容.通常称(8.4.3)为单步法的相容条件.满足相容条件(8.4.3)是可以用单步法求解初值问题(8.1.1)的必要条件.容易验证欧拉法和改进欧拉法均满足相容性条件.一般地，如果单步法有p 阶精度（1p ≥），则其局部截断误差为[]1()()(,(),)()p y x h y x h x y x h O h ϕ++-+=上式两端同除以h ，得()()(,,)()p y x h y x x y h O h hϕ+--=令0h →，如果(,(),)x y x h ϕ连续，则有()(,,0)0y x x y ϕ'-=所以1p ≥的单步法均与问题(8.1.1)相容.由此即得各阶龙格－库塔法与初值问题(8.1.1)相容.定义8.4一种数值方法称为是收敛的，如果对于任意初值0y 及任意固定的(,]x a b ∈，都有lim () ()n h y y x x a nh →==+其中()y x 为初值问题(8.1.1)的精确解.如果我们取消局部化假定，使用某单步法公式，从0x 出发，一步一步地推算到1n x +处的近似值1n y +.若不计各步的舍入误差，而每一步都有局部截断误差，这些局部截断误差的积累就是整体截断误差.定义8.5称111()n n n e y x y +++=-为某数值方法的整体截断误差.其中()y x 为初值问题(8.1.1)的精确解，1n y +为不计舍入误差时用某数值方法从0x 开始，逐步得到的在1n x +处的近似值（不考虑舍入误差的情况下，局部截断误差的积累）.定理8.1设单步法(8.3.14)具有p 阶精度，其增量函数(,,)x y h ϕ关于y 满足利普希茨条件，问题(8.1.1)的初值是精确的，即00()y x y =，则单步法的整体截断误差为111()()p n n n e y x y O h +++=-=证明由已知，(,,)x y h ϕ关于y 满足利普希茨条件，故存在0L >，使得对任意的12,y y 及[,]x a b ∈，00h h <≤，都有1212(,,)(,,)x y h x y h L y y ϕϕ-≤-记1()(,(),)n n n n y y x h x y x h ϕ+=+，因为单步法具有p 阶精度，故存在0M >，使得1111()p n n n R y x y Mh ++++=-≤从而有111111111()()()(,(),)(,,)()(,(),)(,,)n n n n n n n p n n n n n n p n n n n n n e y x y y x y y y Mh y x h x y x h y h x y h Mh y x y h x y x h x y h ϕϕϕϕ+++++++++=-≤-+-≤++--≤+-+-1(1)p nMh hL e +≤++反复递推得11111101110(1)(1)1(1)(1)(1)(1)1(1)p p n n n p n n p n e Mh hL Mh hL e hL hL Mh hL e hL Mh hL e hL+++-+++++⎡⎤≤++++⎣⎦⎡⎤≤+++++++⎣⎦+-≤++因为00()y x y =，即00e =，又(1)n h b a +≤-，于是ln(1)1()(1)(1)b a b a hL n L b a h h hL hL e e --++-+≤+=≤所以()11()p L b a p n M e h e O h L -+⎡⎤≤-=⎣⎦推论设单步法具有p (1p ≥)阶精度，增量函数(,,)x y h ϕ在区域G ：, , 0a x b y h h ≤≤-∞<<+∞≤≤上连续，且关于y 满足利普希茨条件，则单步法是收敛的.当(,)f x y 在区域:,D a x b y ≤≤-∞<<+∞上连续，且关于y 满足利普希茨条件时，改进欧拉法，各阶龙格－库塔法的增量函数(,,)x y h ϕ在区域G 上连续，且关于y 满足利普希茨条件，因而它们都是收敛的.关于单步法收敛的一般结果是：定理8.2设增量函数(,,)x y h ϕ在区域G 上连续，且关于y 满足利普希茨条件，则单步法收敛的充分必要条件是相容性条件(8.4.3).8.4.2稳定性稳定性与收敛性是两个不同的概念，收敛性是在假定每一步计算都准确的前提下，讨论当步长0h →时，方法的整体截断误差是否趋于零的问题.而稳定性则是讨论舍入误差的积累能否对计算结果有严重影响的问题.定义8.6若一种数值方法在节点值n y 上有一个大小为δ的扰动，于以后各节点()m y m n >上产生的偏差均不超过δ，则称该方法是稳定的.我们以欧拉法为例进行讨论.假设由于舍入误差，实际得到的不是n y 而是n n n y y δ=+，其中n δ是误差.由此再计算一步，得到1(,)n n n n y y hf x y +=+把它与不考虑舍入误差的欧拉公式相减，并记111n n n y y δ+++=-，就有[]1(,)(,)1(,)n n n n n n y n nh f x y f x y hf x δδηδ+⎡⎤=+-=+⎣⎦其中y f f y∂=∂.如果满足条件1(,)1y n hf x η+≤，(8.4.4)则从n y 到1n y +的计算，误差是不增的，可以认为计算是稳定的.如果条件(8.4.4)不满足，则每步误差将增大.当0y f >时，显然条件(8.4.4)不可能满足，我们认为问题本身具有先天的不稳定性.当0y f <时，为了满足稳定性要求(8.4.4)，有时h 要很小.一般的，稳定性与方法有关，也与步长h 的大小有关，当然也与方程中的(,)f x y 有关.为简单起见，通常只考虑数值方法用于求解模型方程的稳定性，模型方程为y y λ'=(8.4.5)其中λ为复数.一般的方程可以通过局部线性化转化为模型方程，例如在(,)x y 的邻域(,)(,)(,)()(,)()x y y f x y f x y f x y x x f x y y y '==+-+-+略去高阶项，再作变量替换就得到u u λ'=的形式.对于模型方程(8.4.5)，若Re 0λ>，类似以上分析，可以认为方程是不稳定的.所以我们只考虑Re 0λ<的情形，这时不同的数值方法可能是数值稳定的或者是数值不稳定的.当一个单步法用于试验方程y y λ'=，从n y 计算一步得到1()n n y E h y λ+=(8.4.6)其中()E h λ依赖于所选的方法.因为通过点(,)n n x y 试验方程的解曲线（它满足,()n n y y y x y λ'==）为[]exp ()n n y y x x λ=-，而一个p 阶单步法的局部截断误差在()n n y x y =时有1111()()p n n n T y x y O h ++++=-=，所以有1exp()()()p n n y h E h y O h λλ+-=(8.4.7)这样可以看出()E h λ是h e λ的一个近似值.由(8.4.6)可以看到，若n y 计算中有误差ε，则计算1n y +时将产生误差()E h λε，所以有下面定义.定义8.7如果(8.4.6)式中，()1E h λ<，则称单步法(8.3.14)是绝对稳定的.在复平面上复变量h λ满足()1E h λ<的区域，称为方法(8.3.14)的绝对稳定区域，它与实轴的交称为绝对稳定区间.在上述定义中，规定严格不等式成立，是为了和线性多步法的绝对稳定性定义一致.事实上，()1E h λ=时也可以认为误差不增长.(1)欧拉法的稳定性欧拉法用于模型方程(8.4.5)，得1(1)n n y h y λ+=+，所以有()1E h h λλ=+.所以绝对稳定条件是11h λ+<，它的绝对稳定区域是h λ复平面上以(1,0)-为中心的单位圆，见图8.3.而λ为实数时，绝对稳定区间是(2,0)-.Im()h λRe()h λ2-1-O 图8.3欧拉法的绝对稳定区域(2)梯形公式的稳定性对模型方程，梯形公式的具体表达式为11()2n n n n h y y y y λλ++=++，即11212n nh y y h λλ++=-，所以梯形公式的绝对稳定区域为12112h h λλ+<-.化简得Re()0h λ<，因此梯形公式的绝对稳定区域为h λ平面的左半平面，见图8.4.特别地，当λ为负实数时，对任意的0h >，梯形公式都是稳定的.Im()h λRe()h λO 图8.4梯形公式的绝对稳定区域(3)龙格－库塔法的稳定性与前面的讨论相仿，将龙格－库塔法用于模型方程(8.4.5)，可得二、三、四阶龙格－库塔法的绝对稳定区域分别为211()12h h λλ++<23111()()126h h h λλλ+++<2341111()()()12624h h h h λλλλ++++<当λ为实数时，二、三、四阶显式龙格－库塔法的绝对稳定区域分别为20h λ-<<、2.510h λ-<<、 2.780h λ-<<.例8.5设有初值问题21010101(0)0xy y x x y ⎧'=-≤≤⎪+⎨⎪=⎩用四阶经典龙格－库塔公式求解时，从绝对稳定性考虑，对步长h 有何限制？解对于所给的微分方程有2100,(010)1f x x y xλ∂==-<≤≤∂+在区间[0,10]上，有201010max ||max51t x x λ<<==+由于四阶经典龙格－库塔公式的绝对稳定区间为 2.7850h λ-<<，则步长h 应满足00.557h <<.。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
实验4常微分方程数值解实验目的：1.练习数值积分的计算；2.掌握用MATLAB软件求微分方程初值问题数值解的方法；3.通过实例学习用微分方程模型解决简化的实际问题；4.了解欧拉方法和龙格——库塔方法的基本思想和计算公式，及稳定性等概念。
实验内容：3.小型火箭初始质量为1400kg，其中包括1080kg燃料，火箭竖直向上发射是燃料燃烧率为18kg/s，由此产生32000N的推力，火箭引擎在燃料用尽时关闭。
设火箭上升是空气阻力正比于速度的平方，比例系数为0.4kg/m，求引擎关闭瞬间火箭的高度，速度，加速度，及火箭到达最高点是的高度，速度和加速度，并画出高度，速度，加速度随时间变化的图形。
解答如下：这是一个典型的牛顿第二定律问题，分析火箭受力情况；先规定向上受力为正数建立数学模型：A燃料未燃尽前，在任意时刻（t<60s）火箭受到向上的-F=32000N，向下的重力G=mg，g=9.8，向下的阻力f=kv^2, k=0.4, v表示此时火箭速度；此时火箭收到的合力为F1=（F-mg-f）；火箭的初始质量为1400kg，燃料燃烧率为-18kg/s;此刻火箭质量为m=1400-18*t根据牛顿第二定律知，加速度a=F1/m=(F-mg-f)/(m-r*t)=(32000-(0.4.*v.^2)-9.8.*(1400-18.*t))由此可利用龙格-库塔方法来实现，程序实现如下Function [dx]=rocket[t,x] %建立名为rocket的方程m=1400;k=0.4;r=-18;g=9.8; %给出题目提供的常数值dx=[x(2);(32000-(k*x(2)^2)-g*(m+r*t))/(m+r*t)];%以向量的形式建立方程[a]=(32000-(k*x(2)^2)-g*(m+r*t))/(m+r*t); %给出a的表达式End;ts=0:60; %根据题目给定燃烧率计算出燃料燃尽的时间，确定终点x0=[0,0]; %输入x的初始值[t,x]=ode15s(@rocket,ts,x0); %调用ode15s计算[t,x];h=x(:,1);v=x(:,2);plot(t,x(:,1)),grid; %绘出火箭高度与时间的关系曲线title('h-t');xlabel('t/s');ylabel('h/m'),pause;plot(t,x(:,2)),grid ; %绘出火箭速度与时间的曲线关系title('v-t');xlabel('t/s');ylabel('v/m/s'),pause;a=(32000-(0.4.*v.^2)-9.8.*(1400-18.*t))/(1400-18.*t);plot(t,a),grid; %绘出火箭加速度与时间的曲线关系title('a-t');xlabel('t/s'),ylabel('a/m^2/s'),pause火箭高度随时间变化的曲线火箭速度随时间变化的曲线火箭加速度随时间变化的曲线数据过多，故截取部分如下第一列为时间，第二列为火箭高度，第三列为火箭速度由此可以，在t=60s时，即火箭燃料燃尽瞬间，引擎关闭瞬间，火箭将到达12912m的高度，速度为267,29m，加速度a=0.9m/s^2B燃料燃尽之后，与A 类似，分析受力如下火箭受到向上的F=0向下的重力G=mg，g=9.8，向下的阻力f=kv^2, k=0.4, v表示此时火箭速度；此时火箭收到的合力为F2=（-mg-f）；火箭的初始质量为320kg，恒定根据牛顿第二定律,加速度a=F2/m=-g-0.4v^2/320;程序实现如下function [ dx ] = rocket2( t,x ) %建立以rocket2为名的函数dx=[x(2);-9.8-0.4.*x(2).^2/320]; %以向量的形式建立方程ts=60:120; %给出初始时刻，估计终点时刻x0=[12190,267.26]; %给出x初始值[t,x]=ode15s(@rocket2,ts,x0); %调用ode15s计算[t,x]plot(t,x(:,1)),grid; %绘出火箭高度随时间变化的曲线title('h-t');xlabel('t/s'),ylabel('h/m'),pause;plot(t,x(:,2)),grid; %绘出火箭速度随时间的变化曲线title('v-t');xlabel('t/s'),ylabel('v/m/s'),pause;v=x(:,2);a=-9.8-0.4*v.^2/320; %给出加速度的具体表达式plot(t,a),grid; %绘出火箭加速度随时间变化的曲线title('a-t');xlabel('t/s'),ylabel('a/m^2/s'),pause得到的曲线图形如下火箭高度随时间的变化曲线从图中可以大致看出，最高点在13km左右，火箭速度随时间的变化曲线加速度随时间变化曲线如下数据表格大致如下从图表中可以看出，在71s左右速度到达0，即此时到达最高处，高度为13117m加速度为-9.8m/m/s^2;本题总结：这道题是典型的物理牛顿力学的题目，通过受力的正确分析，可以知道，以[h,v]为向量建立微分方程即可求解，h的微分是速度v，速度v的微分是加速度a解题过程中存在的难点是：取值步长不太容易确定，而且是哪种算法不确定，先用ode15s速度较快，ode23s速度差不太多，其他两种速度较慢，等待时间较长5．一只小船渡过宽为d 的河流，目标是起点A 正对着的另一岸B 点。
已知河水流速为v1 与船在静水中的速度v2之比为k。
（1）建立描述小船航线的数学模型，求其解析解；（2）设d=100m，v1=1m/s，v2=2m/s，用数值解法求渡河所需的时间、任意时刻小船的位置及航行曲线，作图，并与解析解比较。
（3）若流速v1=0，0.5，1.5，2（m/s）, 情况又如何建立数学模型：在任意时刻t，小船位于（x，y），此时速度为v，根据物理中路程与速度的关系，知路程的微分为速度v，由此中，小船在x,y方向上的速度分别为：X：dxdt=v1-v2*√x2+y^2Y:dydt =-v2*2初始条件为dxdt(t0)=v1 *（1）dydt(t0)=v2 *（2）现求其解析解,利用微积分知识*（1）÷*（2）得dxdy =v1∗√x2+y^2−v2y+xy;*(3)将*（3）右端带根号部分分子分母均÷y得dx dy =v1∗√(x/y)2+1−v2+xy;令x/y=p，得到dx/dy=dp/dy*y+p; dx/dy=p;分离变量有y^(-k)=c(√1+p2+p);代入初值可确定当t=0时，y=-d,x=0,p=x/y=0,C=(-d)^(-k)y^(-k)/d^(-k)=√1+x2y^2+x/y;根据隐函数与显函数的关系可得到解析解：x=-d2[(y−d)^(-k)-(yd)^(1+k)]已知d=100m,v1=1m/s,v2=2m/s;先利用龙格库塔方法求解渡河时间，及任意时刻小船的位置（x，y），及航行曲线，与解析解比较此时，k=0.5，d=100用MATLB编写源程序如下：function dx = boat(t,x) %建立名为boat的函数d=-100;v1=1;v2=2; %给定常数值s=(x(1).^2+x(2).^2).^0.5;dx=[v1-v2.*x(1)./s;-v2.*x(2)./s];%用向量形式建立方程ts=0:70; %大致估算，确定终点值，给定步长为”1”x0=[0,-100]; %给出x初始值[t,x]=ode23s(@boat,ts,x0); %调用ode23s计算[t,x];plot(t,x),grid, %绘出x(t),y(t)的函数曲线（图21）gtext('x(t)');gtext('y(t)');pause;plot(t,x(:,1)), %绘出x随时间变化的曲线（图22）grid;xlabel('t/s');ylabel('x/m');pause;plot(t,x(:,2)),grid, %绘出y随时间变化的曲线（图23）title('y-t');xlabel('t/s'),ylabel('y/m');图21图22图23得到数据如下从表格中数据可知，在大约67s时y=0即船到达对岸目的地，为比较，先进行解析解的求解设计程序如下：function x=f(y)k=0.5;x=-0.5.*(-0.01).^k.*y.^(k+1)+0.5.*(-0.01).^(-k).*y.^(-k+1);y=[0:-0.1:-100];for i=0:1:1000;x(:,i+1)=xy(-i/10);endplot(x,y);grid;gtext('x');gtext('y');由此可以看出，由数值解和解析解得到的x-y曲线相差不多，所以可以认为解析解正确改变水流速度v1，只要在原有程序基础上重新复制给v1=0,0.5,1.5,2,同时适当改变终点值即可现实现程序如下A:v1=0，d=-100;v1=0;v2=2;s=(x(1).^2+x(2).^2).^0.5;dx=[v1-v2.*x(1)./s;-v2.*x(2)./s]ts=0:60;x0=[0,-100];[t,x]=ode15s(@boat21,ts,x0);[t,x];plot(x(:,1),x(:,2)),grid,title('y-x');pause;plot(t,x(:,1)),grid;xlabel('t/s');ylabel('x/m');plot(t,x(:,1)),grid;title('x-t');xlabel('t/s'),ylabel('x/m'),pause; plot(t,x(:,2)),grid,title('y-t');xlabel('t/s'),ylabel('y/m');图形如下从此图形中我们可以看到，船并未偏离x=0的点，我们也可以从直观想象中的得到，当水速v1=0时，只要出发时，船头对准目标点，船将一直朝着直线向目的地行进从表格中数据我们也可以很清楚地看到路程与时间是成明显的线性关系的，这是与我们水速为0的必然结果，由此也可以验证我们模型基本正确现改变水速B:v1=0.5程序实现如下d=-100;v1=0.5;v2=2;s=(x(1).^2+x(2).^2).^0.5;dx=[v1-v2.*x(1)./s;-v2.*x(2)./s] ts=0:60;x0=[0,-100];[t,x]=ode15s(@boat22,ts,x0); [t,x];plot(t,x),grid,gtext('x(t)');gtext('y(t)'); pause;plot(t,x(:,1)),grid;xlabel('t/s');ylabel('x/m');plot(t,x(:,1)),grid;title('x-t');xlabel('t/s'),ylabel('x/m'),pause; plot(t,x(:,2)),grid,实现程序过程中发现，终点值设为60，而在53s之后不再有出现，所以可以认为在54s之前就已经达到对岸现在重新改变水速C:v1=1.5d=-100;v1=1.5;v2=2;s=(x(1).^2+x(2).^2).^0.5;dx=[v1-v2.*x(1)./s;-v2.*x(2)./s] ts=0:120;x0=[0,-100];[t,x]=ode15s(@boat23,ts,x0); [t,x];plot(x(:,1),x(:,2)),grid,title('y-x');pause;plot(t,x(:,1)),grid;xlabel('t/s');ylabel('x/m');plot(t,x(:,1)),grid;title('x-t');xlabel('t/s'),ylabel('x/m'),pause; plot(t,x(:,2)),grid,title('y-t');xlabel('t/s'),ylabel('y与v1=0.5s类似，在114s之后不再给出数据，而我们设定的终点值是120，所以可以大致在114s时到达B点现改变水速D:v1=2程序实现如下d=-100;v1=2;v2=2;s=(x(1).^2+x(2).^2).^0.5;dx=[v1-v2.*x(1)./s;-v2.*x(2)./s] ts=0:300;x0=[0,-100];[t,x]=ode15s(@boat24,ts,x0); [t,x];plot(x(:,1),x(:,2)),grid,title('y-x');pause;plot(t,x(:,1)),grid;xlabel('t/s');ylabel('x/m');plot(t,x(:,1)),grid;title('x-t');xlabel('t/s'),ylabel('x/m'),pause; plot(t,x(:,2)),grid,title('y-t');xlabel('t/s'),ylabel('y/m')曲线如下我们可以从此图中看出，当y=0时，x=50，也就是说，船根本没本法到达正对着的B点，而只能到达对岸，我们可以直观地理解假设我们的船在开始时与水速反向，这时船与水的合速度是0，故船无法前进，而根据方程知道，在船尚未到达对岸之前，只有当船的x轴向速度在模式刻超过水速，才有可能克服因水速而在x轴偏离B点的距离，而重新返回，到达B点，而在此给定的速度中，我们可以看到，船速的分量不可能超过水速，因此不可能到达B本题总结：这道题跟课本例题缉私船有较大的相似处，大体是模仿例题而做，所以在画图上显得有点繁琐，如果采用subplot 作图，将会节约一些篇幅，使得作业版面看起来更整洁一点3.两种群的竞争模型如下：其中x(t)，y(t)分别为甲乙两种群的数量；r1，r2为他们的固有增涨率；n1，n2为他们的最大容量。