数学实验第二次作业——常微分方程数值求解
常微分方程数值解实验
有大量的常微分方程虽然从理论上讲，其解是存在的，但我们却无 法求出其解析解，此时，我们需要寻求方程的数值解，在求常微分方程 数值解方面，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))得
常微分方程数值解实验报告
常微分方程数值解实验报告学院：数学与信息科学专业：信息与计算科学姓名：郑思义学号：201216524课程：常微分方程数值解实验一：常微分方程的数值解法1、分别用Euler 法、改进的Euler 法（预报校正格式）和S —K 法求解初值问题。
（h=0.1）并与真解作比较。
⎩⎨⎧=++-=10(1y'）y x y 1.1实验代码：%欧拉法function [x,y]=naeuler(dyfun,xspan,y0,h)%dyfun 是常微分方程，xspan 是x 的取值范围，y0是初值，h 是步长 x=xspan(1):h:xspan(2); y(1)=y0;for n=1:length(x)-1y(n+1)=y(n)+h*feval(dyfun,x(n),y(n)); end%改进的欧拉法function [x,m,y]=naeuler2(dyfun,xspan,y0,h)%dyfun 是常微分方程，xspan 是x 的取值范围，y0是初值，h 是步长。
%返回值x 为x 取值，m 为预报解，y 为校正解 x=xspan(1):h:xspan(2); y(1)=y0;m=zeros(length(x)-1,1); for n=1:length(x)-1 k1=feval(dyfun,x(n),y(n)); y(n+1)=y(n)+h*k1; m(n)=y(n+1);k2=feval(dyfun,x(n+1),y(n+1)); y(n+1)=y(n)+h*(k1+k2)/2; end%四阶S —K 法function [x,y]=rk(dyfun,xspan,y0,h)%dyfun 是常微分方程，xspan 是x 的取值范围，y0是初值，h 是步长。
x=xspan(1):h:xspan(2); y(1)=y0;for n=1:length(x)-1 k1=feval(dyfun,x(n),y(n));k2=feval(dyfun,x(n)+h/2,y(n)+(h*k1)/2); k3=feval(dyfun,x(n)+h/2,y(n)+(h*k2)/2); k4=feval(dyfun,x(n)+h,y(n)+h*k3);y(n+1)=y(n)+(h/6)*(k1+2*k2+2*k3+k4);end%主程序x=[0:0.1:1];y=exp(-x)+x;dyfun=inline('-y+x+1');[x1,y1]=naeuler(dyfun,[0,1],1,0.1);[x2,m,y2]=naeuler2(dyfun,[0,1],1,0.1);[x3,y3]=rk(dyfun,[0,1],1,0.1);plot(x,y,'r',x1,y1,'+',x2,y2,'*',x3,y3,'o');xlabel('x');ylabel('y');legend('y为真解','y1为欧拉解','y2为改进欧拉解','y3为S—K解','Location','NorthWest');1.2实验结果：x 真解y 欧拉解y1 预报值m 校正值y2 S—K解y30.0 1.0000 1.0000 1.0000 1.00000.1 1.0048 1.0000 1.0000 1.0050 1.00480.2 1.0187 1.0100 1.0145 1.0190 1.01870.3 1.0408 1.0290 1.0371 1.0412 1.04080.4 1.0703 1.0561 1.0671 1.0708 1.07030.5 1.1065 1.0905 1.1037 1.1071 1.10650.6 1.1488 1.1314 1.1464 1.1494 1.14880.7 1.1966 1.1783 1.1945 1.1972 1.19660.8 1.2493 1.2305 1.2475 1.2500 1.24930.9 1.3066 1.2874 1.3050 1.3072 1.30661.0 1.3679 1.3487 1.3665 1.3685 1.36792、选取一种理论上收敛但是不稳定的算法对问题1进行计算，并与真解作比较。
实验报告——常微分方程的数值解法
实验报告实验项目名称常微分方程的数值解法实验室数学实验室所属课程名称微分方程数值解实验类型上机实验实验日期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 法的精度更高。
常微分方程数值解实验
多步法，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日
数学实验——常微分方程数值解
实验4 常微分方程数值解分1 黄浩 43一、实验目的1.掌握用MATLAB软件求微分方程初值问题数值解的方法；2.通过实例学习用微分方程模型解决简化的实际问题；3.了解欧拉方法和龙格-库塔方法的基本思想和计算公式，及稳定性等概念。
二、实验内容1.《数学实验》第一版（问题2）问题叙述：小型火箭初始重量为1400kg，其中包括1080kg燃料。
火箭竖直向上发射时燃料燃烧率为18kg/s，由此产生32000N的推力，火箭引擎在燃烧用尽时关闭。
设火箭上升时空气阻力正比于速度的平方，比例系数为0.4kg/m，求引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点时的高度和加速度，并画出高度、速度、加速度随时间变化的图形。
模型转换及实验过程：（一）从发射到引擎关闭设火箭总质量为m，上升高度为h，瞬时速度为v，瞬时加速度为a，由燃料燃烧时间t=60s，可列如下的方程组：t∈[0,60]−9.8因此，上述方程为二元常微分方程组，选择t为自变量，h和v为因变量进行分析。
初值条件：h(0)=0 ,v(0)=0对上述模型，使用ode45()函数求数值解（程序见四.1、四.2），结果如下：由上表可知，引擎关闭瞬间，火箭的高度为12189.78m，速度为267.26m/s，加速度为0.9170m/s2，火箭至此已飞行60s而高度、速度、加速度随时间的变化曲线如下：(二)从引擎关闭到最高点设引擎关闭时，t1=0，由上一问的结果可知，h(0)=12189.78m，v(0)= 267.26m/s，m=320kg，则可列二元常微分方程组如下：9.8因此，可选择t1为自变量，h、v为因变量进行分析（程序见四.3、四.4），实验结果如下：由上表可知，当t1∈[11,12]时，v(t1)有零点，即该区间内某时刻火箭达到最高点。
再进行更细致的实验（程序略），设步长为0.01，观察该区间内v(t1)的零点，如下表所示：可以看出，当t1=11.30s，即总时间t=71.30s时，火箭达到最高点，高度为13115.36m，加速度为-9.8m/s2。
常微分方程的数值求解
常微分方程的数值求解在数学中，常微分方程是一类重要的数学模型，通常用来描述物理、化学、生物等自然现象中的变化规律。
对于一些复杂的微分方程，无法通过解析方法进行求解，这时候就需要借助数值方法来近似求解。
本文将介绍常微分方程的数值求解方法及其应用。
一、数值求解方法常微分方程的数值求解方法主要包括欧拉法、改进的欧拉法、龙格-库塔法等。
欧拉法是最简单也是最常用的数值求解方法，其基本思想是根据微分方程的导数近似求解下一个时间步上的解，并通过逐步迭代来得到整个解的数值近似。
改进的欧拉法在欧拉法的基础上做出了一定的修正，提高了数值求解的精度。
而龙格-库塔法则是一种更加精确的数值求解方法，通过考虑多个点的斜率来进行求解，从而减小误差。
二、应用领域常微分方程的数值求解方法在科学研究和工程实践中有着广泛的应用。
在物理学中，通过数值求解微分方程可以模拟天体运动、粒子运动等现象；在生物学领域，可以模拟生物种群的增长和变化规律；在工程领域，可以通过数值求解微分方程来设计控制系统、优化结构等。
三、实例分析以一个简单的一阶常微分方程为例：dy/dx = -y，初始条件为y(0) = 1。
我们可以用欧拉法来进行数值求解。
将时间间隔取为0.1，通过迭代计算可以得到y(1)的近似值为0.367。
而利用改进的欧拉法或者龙格-库塔法可以得到更加精确的数值近似。
这个例子展示了数值方法在解决微分方程问题上的有效性。
四、总结常微分方程是求解自然界中变化规律的重要数学工具，而数值方法则是解决一些难以解析求解的微分方程的有效途径。
通过本文的介绍，读者可以了解常微分方程的数值求解方法及其应用，希望可以对相关领域的研究和实践有所帮助。
至此，关于常微分方程的数值求解的文章正文部分结束。
数学实验第二次作业——常微分方程数值求解
实验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 点。
常微分方程数值解法
i.常微分方程初值问题数值解法i.1 常微分方程差分法考虑常微分方程初值问题：求函数()u t 满足(,), 0du f t u t T dt=<≤ （i.1a ） 0(0)u u = (i.1b)其中(,)f t u 是定义在区域G : 0t T ≤≤, u <∞上的函数，0u 和T 是给定的常数。
我们假设(,)f t u 对u 满足Lipschitz 条件，即存在常数L 使得121212(,)(,), [0,]; ,(,)f t u f t u L u u t T u u -≤-∀∈∈-∞∞ （i.2） 这一条件保证了（i.1）的解是适定的，即存在，唯一，而且连续依赖于初值0u 。
通常情况下，（i.1）的精确解不可能用简单的解析表达式给出，只能求近似解。
本章讨论常微分方程最常用的近似数值解法--差分方法。
先来讨论最简单的Euler 法。
为此，首先将求解区域[0,]T 离散化为若干个离散点：0110N N t t t t T -=<<<<= （i.3） 其中n t hn =，0h >称为步长。
在微积分课程中我们熟知，微商（即导数）是差商的极限。
反过来，差商就是微商的近似。
在0t t =处，在（i.1a ）中用向前差商10()()u t u t h -代替微商du dt ，便得 10000()()(,())u t u t hf t u t ε=++如果忽略误差项0ε，再换个记号，用i u 代替()i u t 便得到1000(,)u u hf t u -=一般地，我们有1Euler (,), 0,1,,1n n n n u u hf t u n N +=+=-方法： （i.4） 从(i.1b) 给出的初始值0u 出发，由上式可以依次算出1,,N t t 上的差分解1,,N u u 。
下面我们用数值积分法重新导出 Euler 法以及其它几种方法。
为此，在区间1[,]n n t t +上积分常微分方程（i.1a ），得11()()(,())n n t n n t u t u t f t u t dt ++=+⎰ （i.5）用各种数值积分公式计算（i.5）中的积分，便导致各种不同的差分法。
数值分析实验报告之常微分方程数值解
数学与计算科学学院实验报告
实验项目名称常微分方程数值解
所属课程名称数值方法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. 根据实际要求确定要研究的量(自变量、未知函数、必要的参数等)并确定坐标系。
2. 找出这些量所满足的基本规律(物理的、几何的、化学的或生物学的等等)。
3. 运用这些规律列出方程和定解条件。
基本模型1. 发射卫星为什么用三级火箭2. 人口模型3. 战争模型4. 放射性废料的处理通常需要求出方程的解来说明实际现象，并加以检验。
如果能得到解析形式的解固然是便于分析和应用的，但是我们知道，只有线性常系数微分方程，并且自由项是某些特殊类型的函数时，才可以得到这样的解，而绝大多数变系数方程、非线性方程都是所谓“解不出来” 的于是对于用微分方程解决实际问题来说，数值解法就是一个十分重要的手段。
1. 改进Euler 法：2. 龙格—库塔( Runge—Kutta )方法：【源程序】1. 改进Euler 法：function [x,y]=eulerpro(fun,x0,x1,y0,n);%fun 为函数，(xO, x1)为x 区间，yO 为初始值，n 为子区间个数if nargin<5,n=5O;endh=(x1-xO)/n;x(1)=xO;y(1)=yO;for i=1:nx(i+1)=x(i)+h;y1=y(i)+h*feval(fun,x(i),y(i));y2=y(i)+h*feval(fun,x(i+1),y1);y(i+1)=(y1+y2)/2;end调用command 窗口f=i nlin e('-2*y+2*x A2+2*x') [x,y]=eulerpro(f,O,,1,1O)2x +2x , (0 < x < , y(0) = 1求解函数y'=-2y+22. 龙格—库塔( Runge—Kutta )方法：[t,y]=solver('F',tspan ，y0)这里solver为ode45, ode23, ode113，输入参数F是用M文件定义的微分方程y'= f (x, y)右端的函数。
实验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); }}。
数学建模-- 常微分方程数值解及实验
注意:
1、在解n个未知函数的方程组时，x0和x均为n维向量， m-文件中的待解方程组应以x的分量形式写成. 2、使用Matlab软件求数值解时，高阶微分方程必须 等价地变换成一阶微分方程组.
例
d 2x dx 2 x 0 2 ( x 1) dt dt x (0 ) 2; x '(0 ) 0
实际应用时，与欧拉公式结合使用：
0 y i( 1) y i hf ( x i , y i ) h ( k 1) (k ) y i 1 y i [ f ( x i , y i ) f ( x i 1 , y i 1 )] k 0 ,1, 2 , 2
•欧拉法是一阶公式，改进的欧拉法是二阶公式。
•龙格-库塔法有二阶公式和四阶公式。 •线性多步法有四阶阿达姆斯外插公方程的数值解
[t，x]=solver（’f’,ts,x0,options）
自变 量值 函数 值
ode45 ode23 ode113 ode15s ode23s
1、建立m-文件vdp1.m如下： function dy=vdp (t,y) dy=zeros(2,1); dy(1)=y(2); dy(2)=(1-y(1)^2)*y(2)-y(1); 2、取t0=0，tf=20，输入命令： [T,Y]=ode15s('vdp1',[0 20],[2 0]); plot(T,Y(:,1),'-')
y ( x i 1 ) y ( x i )
xi 1
f ( t , y ( t )) dt
x i 1 x i 2
xi
[ f ( x i , y ( x i )) f ( x i 1 , y ( x i 1 ))]
常微分方程数值解---实验报告
琼州学院实验报告课程名称：__________ ___开课学期：____ ______院（部）:__理工学院_________开课实验室：________ __学生姓名:__梁小叶_______ __专业班级:________ __学号:__________常微分方程数值解---实验报告一 实验目的： 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
常微分方程的数值解法
常微分方程的数值解法常微分方程的数值解法，这个话题一听就让人感觉有点儿“高大上”，是不是？别急，咱们今天就把这个看似复杂的东西聊得轻松点。
想象一下，你的生活中随时都在用到数学，走路时的速度、车子的加速，甚至你做的美食，背后其实都藏着一些微分方程的身影。
说白了，常微分方程就是描述一个变量随另一个变量变化的关系。
听起来是不是有点儿抽象？咱们可以用一些简单的例子来说明。
比如说，想象一下你在家里的花园里浇水，植物的生长就可以用微分方程来描述。
你浇水越多，植物长得越快，没错吧？这时候，就有个问题来了：你浇水的频率和量应该是多少呢？这就涉及到微分方程了。
可是，现实生活中，咱们并不能总是拿个公式就能算出结果。
假如你想要把这个过程模拟出来，找出植物生长的规律，这时候数值解法就派上用场了。
数值解法听起来很复杂，但实际上它的核心思想就是通过一些巧妙的方法把这些难以直接求解的微分方程转化为我们能够计算的形式。
就像是你要吃一个巨大的蛋糕，如果直接咬下去，肯定会噎着；所以，聪明的你会选择先把蛋糕切成小块，再慢慢享用。
这就是数值解法的思路。
它把复杂的事情分解成简单的小步骤，让我们能一个一个地解决。
常见的数值解法有很多，比如欧拉法、龙格库塔法。
这些名字听起来像是外星人的语言，其实很简单。
以欧拉法为例，就像在慢慢走路，每一步都非常稳当。
你先确定一个初始值，然后根据微分方程的斜率来推算下一步的值。
一步一步走，虽然速度慢点儿，但稳妥得很，不容易出错。
只要你把每一步都计算好，最终就能到达目标。
再说说龙格库塔法，这个名字听起来像是在做科学实验，其实它是一种更高效的方法。
就好比你打游戏，普通模式下你需要慢慢磨练，而龙格库塔法就像是开了个外挂，让你更快地达到目标。
这个方法通过在每一步中做多次的预测和修正，来提高精度，减少误差。
简单来说，就是“多看看再决定”，这样能让你走得更稳、更远。
咱们在应用这些方法时，当然得注意一些细节。
比如说，步长的选择。
常微分方程数值解法实验报告
常微分方程数值解法实验报告
一、实验目的
1、熟悉各种初值问题的算法，编出算法程序；
2、明确各种算法的精度与所选步长有密切关系；通过计算更加了解各种算法的优越性。
二、实验题目
1、根据初值问题数值算法，分别选择二个初值问题编程计算；
2、试分别取不同步长，考察某节点jx处数值解的误差变化情况；
3、试用不同算法求解某初值问题，结果有何异常；
4、分析各个算法的优缺点。
三、实验原理与理论基础
（一）欧拉法算法设计对常微分方程初始问题(6-1)(6-2)用数值方法求解时，我们总是认为(6-1)、(6-2)的解存在且唯一。
欧拉法是解初值问题的最简单的数值方法。
从(6-2)式由于y(x 0)=y 0已给定，因而可以算出),()('0 0 0y x f x y）=。
记 ,n , , i x y yi i1 0 ) ( = =从而我们可以取。
数学实验报告——利用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 自带龙格-库塔公式的结果稍近。
常微分方程数值解法
int i; x=x0; y=y0; h=(xn-x0)/n; for(i=1;i<=n;i++) { y=y+h*func(x,y); x=x0+i*h; } return(y); } void main() { float x0,xn,y0,e; int n; printf("\ninput n:\n"); scanf("%d",&n); printf("input x0,xn:\n"); scanf("%f %f",&x0,&xn); printf("input y0:\n"); scanf("%f",&y0); e=euler(x0,xn,y0,n); printf("y(%f)=%6.6f",y0,e); } 运行结果为：
实习、改进的欧拉法的基本思想和过程； （3） 实习前复习欧拉法、改进的欧拉法的计算步骤。
实习设备：
（1） 硬件设备：单机或网络环境下的微型计算机一台; （2） 软件设备：DOS3.3以上操作系统，Turbo C2.0编译器。
实习内容：
欧拉法与改进欧拉法 （1） 分别使用欧拉法与改进欧拉法求解初值问题
的数值解。 （2） 要求： 请写出分别使用欧拉法与改进欧拉法程序的运行结果。 欧拉法代码： #include "stdio.h" #include "conio.h" float func(float x,float y) { return(-2*x*y); } float euler(float x0,float xn,float y0,int n) { float x,y,h;
- 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为他们的最大容量。