数学实验第二次作业任务常微分方程数值求解
常微分方程数值解实验
有大量的常微分方程虽然从理论上讲，其解是存在的，但我们却无 法求出其解析解，此时，我们需要寻求方程的数值解，在求常微分方程 数值解方面，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 <<.。
计算方法实验报告-常微分方程的数值解法
15.1000 50.1224 50.1224
15.2000 50.1332 50.1332
15.3000 50.1436 50.1436
15.4000 50.1535 50.1535
disp([T',V',df']);
%画图观察效果
figure;
plot(T,df,'k*',T,V,'--r');
grid on;
title('速度v与时间t的函数曲线');
legend('准确值','四阶经典R-K法');
结果如下所示：
步长四阶经典R-K法准确值
0 0 0
0.1000 0.9799 0.9799
5.9000 41.1708 41.1708
6.0000 41.4919 41.4919
6.1000 41.8028 41.8028
6.2000 42.1039 42.1039
6.3000 42.3954 42.3954
6.4000 42.6775 42.6775
6.5000 42.9504 42.9504
11.5000 49.2667 49.2667
11.6000 49.3097 49.3097
11.7000 49.3511 49.3511
11.8000 49.3909 49.3909
11.9000 49.4292 49.4292
12.0000 49.4661 49.4661
12.1000 49.5016 49.5016
（2）某跳伞者在t=0时刻从飞机上跳出，假设初始时刻的垂直速度为0，且跳伞者垂直下落。已知空气阻力为F=cv2，其中c为常数，v为垂直速度，向下方方向为正。写出此跳伞者的速度满足的微分方程；若此跳伞者的质量为M=70kg，且已知c=0.27kg/m，利用四阶龙格-库塔公式计算t<=20s的速度（取h=0.1s）
数值分析实验报告之常微分方程数值解
数学与计算科学学院实验报告
实验项目名称常微分方程数值解
所属课程名称数值方法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．指导教师评语及成绩：指导教师依据学生的实际报告内容，给出本次实验报告的评价。
数值分析中的常微分方程数值求解
数值分析中的常微分方程数值求解常微分方程是自然科学中一类最为普遍的数学模型，涉及到热力学、物理、化工等多个领域。
然而，解常微分方程并非易事。
尤其是当我们面对一些复杂、非线性、多维的方程组时，常微分方程数值求解成为了一个十分关键的问题。
因此，数值求解方法成为了常微分方程研究中的重要组成部分。
本文将介绍一些数值解常微分方程的常见方法和应用。
1. 一般线性方法一般线性方法（general linear methods）是经典的常微分方程数值解法之一。
它以一种特殊的形式给出步进公式：$$ y_{n+1}=\sum_{i=0}^{s-1}\alpha_i y_{n-i}+h\sum_{i=0}^{s-1}\beta_i f(t_{n-i},y_{n-i}) $$ 其中，$y_{n}$为第$n$步的项值，$f(t_n,y_n)$为时间$t_n$处函数$y(t)$的导数。
$\alpha_i$和$\beta_i$是常数，可以通过确定如下特征方程来选择：$$ \sum_{i=0}^{s-1}\alpha_i\ lambda^{i}=0,~(\lambda\in C) $$ 与此同时，也可以通过选择$\beta_i$来使方法达到一定的准确性和稳定性。
2. Runge-Kutta方法比一般线性方法更为流行的方法是Runge-Kutta方法。
通常附加一个或多个修正以获得更好的数值稳定性和误差控制。
第1阶Runge-Kutta方法仅使用导数$f(t_n,y_n)$估算下一个项的值：$$y_{n+1}=y_n+hf(t_n,y_n)$$ 许多高阶方法可以使用中间的“插值”来更准确地估计下一个步骤：$$y_{n+1}=y_n+h\sum_{i=1}^kb_ik_i$$$$k_i=f(t_n+c_ih,y_n+h\sum _{j=1}^{i-1}a_{ij}k_j)$$ $k_i$是第$i$台车的估计值，$a_{ij}$和$b_i$在经典Runge-Kutta方法和其他变体中具有不同的取值。
实验2常微分方程数值解法
实验2 常微分方程数值解法实验目的与要求1）熟悉求解常微分方程初值问题的有关方法和理论，主要是改进 欧拉法和四阶龙格库塔法2）会编制改进欧拉法的计算程序；3）针对实验内容编制的程序能够正确调试并运行所需结果实验内容：请用改进欧拉法和四阶龙格库塔算法解算法介绍 1、改进欧拉算法概要 20011112111205(0)2(,)()(,)[(,)(,)](,)[(,)(b an i i i i hi i i i i i i i i i h i i i i i y xy x y y f x y a x b y x y y y hf x y y y f x y f x y y y hf x y y y f x y f x 解一阶常微分方程初值问题将区间[a,b]作n 等分，取步长h=欧拉公式为梯形公式为改进欧拉法的公式为-+++++++ì¢ï=-＃ïíï=ïîì¢=＃ïïíï=ïî=+=++=+=++111,)](,)(,)1()2i p i i i c i i i i p c y y y hf x y y y hf x y y y y 或表示为+++ìïïíïïîìïïï=+ïïï=+íïïïï=+ïïî2、龙格库塔算法概要算法函数示例：void ModEuler(float(*f)(float,float),float x0,float y0,float xn,int n){int i;float yp,yc,x=x0,y=y0,h=(xn-x0)/n;printf("x[0]=%f\ty[0]=%f\n",x,y);for(i=1;i<=n;i++){yp=y+h*(*f)(x,y);x=x0+i*h;yc=y+h*(*f)(x,yp);y=(yp+yc)/2;printf("x[%d]=%f\ty[%d]=%f\n",i,x,i,y);}}void Runge_Kutta(float(*f)(float x,float y),float a,float b,float y0,int n){float x=a,y=y0,K1,K2,K3,K4;float h=(b-a)/N;int i;printf("x[0]=%f\ty[0]=%f\n",x,y);for(i=1;i<=n;i++){K1=(*f)(x,y);K2=(*f)(x+h/2,y+h*K1/2);K3=(*f)(x+h/2,y+h*K2/2);K4=(*f)(x+h,y+h*K3);上各个节点初的近似值，在区间出发，可得未知函数，由初值取步长为库塔公式常用的四阶龙格对于初值问题：][a )(),()2,2()2,2(),()22(6)(...........).........,(034231214321100'b x y y h hk y h x f k k h y h x f k k h y h x f k y x f k k k k k h y y y x y b x a y x f y i i i i i i i i i i ⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧++=++=++==++++=-⎩⎨⎧=≤≤=+y=y+h*(K1+2*K2+2K3+K4)/6;x=a+i*h;printf("x[%d]=%f\ty[%d]=%f\n",i,x,i,y); }}。
常微分方程的数值解法实验报告
常微分方程的数值解法实验报告实验报告：常微分方程的数值解法摘要：常微分方程（ODE）是描述动力学系统中物理量随时间变化的数学方程，广泛应用于自然科学和工程领域。
然而，对于一些复杂的非线性ODE，很难找到解析解。
因此，我们需要数值解法来求解这些方程。
本实验报告将介绍四种常见的常微分方程数值解法，分别是欧拉法、改进的欧拉法、四阶龙格-库塔法和自适应步长的龙格-库塔法，并通过数值实验比较它们的精度和效率。
1.引言在实际问题中，许多物理量的变化规律可以由常微分方程描述。
然而，对于复杂的非线性ODE，很难找到解析解。
因此，为了解决这类问题，我们需要借助数值方法来求解。
2.方法本实验采用四种常见的常微分方程数值解法：欧拉法、改进的欧拉法、四阶龙格-库塔法和自适应步长的龙格-库塔法。
（1）欧拉法是最简单的数值解法，通过将微分方程转化为差分方程，使用离散的步长来近似微分方程。
（2）改进的欧拉法在欧拉法的基础上进行了改进，使用预估-校正的方法来提高精度。
（3）四阶龙格-库塔法是一种经典的数值解法，通过利用不同步长处的斜率来近似微分方程，具有较高的精度。
（4）自适应步长的龙格-库塔法是在四阶龙格-库塔法的基础上改进而来的，根据步长的大小自适应地选择不同的步长，同时保证精度和效率。
3.实验设计为了比较这四种数值解法的精度和效率，我们设计了两个实验。
实验一是求解一阶常微分方程：dy/dx = -2x，初始条件y(0) = 1，解析解为y = 1 - x^2、实验二是求解二阶常微分方程：d^2y/dx^2 + y = 0，初始条件y(0) = 0，dy/dx(0) = 1，解析解为y = sin(x)。
4.结果与分析实验一中，比较四种数值解法在不同步长下的近似解和解析解，计算其误差。
实验结果表明，四阶龙格-库塔法和自适应步长的龙格-库塔法具有较高的精度，而欧拉法和改进的欧拉法的精度较低。
实验二中，我们比较四种数值解法在不同步长下的近似解和解析解，并计算其误差。
常微分方程数值解---实验报告
琼州学院实验报告课程名称：__________ ___开课学期：____ ______院（部）:__理工学院_________开课实验室：________ __学生姓名:__梁小叶_______ __专业班级:________ __学号:__________常微分方程数值解---实验报告一 实验目的： 1．掌握用MATLAB 求微分方程初值问题数值解的方法；2．通过实例学习微分方程模型解决简化的实际问题；3．了解欧拉方法和龙格库塔方法的基本思想。
二 实验内容:用欧拉方法和龙格库塔方法求下列微分方程初值问题的数值解，画出解的图形，对结果进行分析比较22'2,(1)(01),322;(0)1',(2)(010).(0)0(0)1x y y x x y e x y y x y x y y =+⎧≤≤=--⎨=⎩⎧=-≤≤⎨==⎩精确解或三 问题分析：怎样设计程序？流程图？变量说明？能否将某算法设计成具有形式参数的函数形式？【程序如下】:function f=f(x,y)f=y+2*x;clc;clear;a=0;b=1; %求解区间[x1,y_r]=ode45('f',[a b],1); %调用龙格库塔求解函数求解数值解；%% 以下利用Euler 方法求解y(1)=1;N=100;h=(b-a)/N;x=a:h:b;for i=1:Ny(i+1)=y(i)+h*f(x(i),y(i));endfigure(1)plot(x1,y_r,'r*',x,y,'b+',x,3*exp(x)-2*x-2,'k-');%数值解与真解图title('数值解与真解图');legend('RK4','Euler','真解');xlabel('x');ylabel('y');figure(2)plot(x1,abs(y_r-(3*exp(x1)-2*x1-2)),'k-');%龙格库塔方法的误差title('龙格库塔方法的误差')xlabel('x');ylabel('Error');figure(3)plot(x,abs(y-(3*exp(x)-2*x-2)),'r-')%Euler 方法的误差title('Euler 方法的误差')xlabel('x');ylabel('Error');【运行结果如下】：5总结体会：自己写了哦。
数学实验第2次作业_常微分方程数值解
22
2804.62262791368
214.164243500080
3.80623366610031
23
3020.62050869921
217.768226726238
3.42270756122453
24
3240.03823406982
221.009658622529
7.65487748743591
16
1604.92980511976
182.532445472022
6.99682845896217
17
1790.85149669532
189.206450509489
6.36706376524626
18
1983.13803757204
195.266573002515
1.09643186923185
40
7031.88138581493
248.037936243021
1.06967218025673
41
7280.45011340898
249.095404105782
1.04741271358073
42
7530.06554547068
250.132761670182
75
13073.2000008109
29.4639811853995
-10.9065
76
13088.9460668350
23.0648815008400
-10.4835
77
13100.9213897508
16.8837287250651
-10.1729
78
13109.2374618597
数学实验报告——利用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 自带龙格-库塔公式的结果稍近。
数学实验第二次作业任务常微分方程数值求解
数学实验第⼆次作业任务常微分⽅程数值求解实验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 点。
- 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为他们的最大容量。