数学实验课程设计常微分方程数值解
常微分方程数值解实验
有大量的常微分方程虽然从理论上讲，其解是存在的，但我们却无 法求出其解析解，此时，我们需要寻求方程的数值解，在求常微分方程 数值解方面，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 法的精度更高。
常微分方程数值解及实验
改进的欧拉法
实际应用时，结合欧拉公式采用迭代法提高精度：
0 yi(1) yi hf ( xi , yi ) h ( k 1) k yi 1 yi [ f ( xi , yi ) f ( xi 1 , yi(1) 当 yi(11) yi(1) 时， yi 1 yi k11） 取 。
y1 ' y2 2 y2 ' (1 y1 ) y2 y1 y (0) 2, y (0) 0 2 1
1、建立m-文件vdp.m如下： function dy=vdp(t,y) dy=zeros(2,1); dy(1)=y(2); dy(2)=(1-y(1)^2)*y(2)-y(1);
（二）建立数值解法的一些途径
设 xi 1 xi h, i 0,1, 2, n 1, 可用以下离散化方法求解微分方程： y ' f ( x, y ) y ( x0 ) y0
1.欧拉方法 若步长h较小，则有 向前欧拉法 向后欧拉法
yi 1 yi hf ( xi , yi ) i 0,1, 2, , n -1 y0 y ( x0 )
yi 1 yi hf ( xi 1 , yi 1 ) i 0,1, 2, , n -1 y0 y ( x0 )
梯形法
f ( xi , yi ) f ( xi 1 , yi 1 ) yi 1 yi h i 0,1, 2, , n -1 2 y0 y ( x0 )
注意: 1.在解n个未知函数的方程组时，x0和x均为n维向量， m-文件中的待解方程组应以x的分量形式写成.
% f.m function x=f(t) .......
常微分方程数值解实验
多步法，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。
实验七、常微分方程数值解法
实验七 常微分方程数值解法实验目的1、 掌握常微分方程数值解法中欧拉法和改进欧拉法的基本原理和算法；2、 培养Matlab 数学软件的编程与上机调试能力.实验课时4课时实验准备 初值问题()(),,dyf x y a x bdxy a η⎧=≤≤⎪⎨⎪=⎩的解析解法在常微分方程中有介绍,但是实际中遇到的常微分方程一般很难得到解析解,所以我们一般用数值方法求出它的数值解. 数值解法的结果是自变量x 在一系列离散点,,21n x x x <<处的函数()x f 的近似值 ,,,,21n y y y ,一般情况下取步长n n x x h -=+1为定数.1.欧拉法利用一阶微分公式()()1nn n x x y x y x dy dxh+=-≈可得差分方程()1,,0,1,...,1n n n n y y hf x y n N +=+=-其中, 0y η=, ()n n y y x ≈,则得到初值问题的数值解.2.后退欧拉法利用一阶微分公式()()11n n n x x y x y x dy dxh++=-≈可得差分方程()1110,,0,1,...,1n n n n y y hf x y n N y η+++=+=-⎧⎪⎨=⎪⎩ 此为后退欧拉法.4. 梯形求积公式利用数值积分()()()()()()()1111,,,2n nx n n n n n n x h y x y x f x y dx f x y x f x y x ++++⎡⎤-=≈+⎣⎦⎰可得梯形公式()()1110,,,0,1,...,12n n n n n n h y y f x y f x y n N y η+++⎧=++=-⎡⎤⎪⎣⎦⎨⎪=⎩5. 改进欧拉公式一般来讲,隐式格式要比显示格式具有较好的数值稳定性,所以常常被采用,但是为了避免隐式格式的计算量,一般采用预估-校正的方法.对于梯形公式,用欧拉公式预估再用梯形公式校正,可得改进欧拉公式:()()()()01111,,,2n n n n n n n n n n y y hf x y hy y f x y f x y ++++⎧=+⎪⎨=++⎡⎤⎪⎣⎦⎩实验算例1. 欧拉法利用matlab 中的一个循环语句即可实现欧拉法中的利用第一个节点计算随后所有节点例1. 用欧拉法求解初值问题, 并画出解图形.()1,0101y y x x y '=-++≤≤⎧⎪⎨=⎪⎩ >>输出图像为另外,还可以改变输入的参数n,加密节点,得到更精确的数值解.对于例1的结果为输出图像为欧拉法和改进欧拉法的效果比较图为实验习题1. 利用欧拉法程序,求解初值问题（注意：若理论解未知，解曲线画不出来，应删掉相关程序语句）(1)()2sin ,1211y y y x x y '⎧=--≤≤⎪⎨=⎪⎩ (2)(),0101yy xe x y '⎧=≤≤⎪⎨=⎪⎩ 并将欧拉法求出的折线和改进欧拉法求出的折线放在同一个图中，比较效果。
数学实验课程设计常微分方程数值解
数学实验报告1.题目：某容器盛满水后，底端直径为d0的小孔开启（如图1），根据水力学知识，当水面高度为h时，谁从小孔中流出的速度为v=0.6*（g*h）^0.5（其中g为重力加速度，0.6问哦小孔收缩系数）1）若容器为倒圆锥形（如图1），现测得容器高和上底直径都为1.2m，小孔直径d为3cm，为水从小孔中流完需要多少时间；2min时水面高度是多少。
2）若容器为倒葫芦形（如图2），现测得容器高1.2m，小孔直径d为3cm，由底端（记x=0）向上每隔0.1m测出容器的直径D（m）如表1所示，问水从小孔中流完需要多少时间；2min时水面高度是多少。
X/m00.10.20.30.40.50.60.70.80.9 1.0 1.1 1.2D/m0.030.050.080.140.190.330.450.680.981.11.21.131.02.分析：由题知，水从小孔中流出，不仅与容器有关，还与水流速度v=0.6*（2*g*h）^0.5有关。
第一小题容器是圆锥形，比较规则，但是由于水不断从小孔流出，容器中水的高度是不断变化的，水流速度没有一定的公式，所以要用到微积分解决。
由（1）知，水面直径等于水深。
水深为h时，流量为0.6(π/4)d^2*(gh)^0.5，0.6*(g*h)^(0.5)*π*(d0/2)^2*dt=π/4*h^2*dh则水深下降dh 所需时间 ：dt=-[(π/4)h^2*dh]/[0.6(π/4)d^2*(gh)^0.5]=-[h^1.5*dh]/[0.6d^2*(g)^0.5]水深由1.2m 至0定积分得水从小孔流完的时间：T （其中已知d=0.03m ，g=9.8m*s(-2）对于第二问：设两分钟（120S ）后水深为X m ，由dt=-[(π/4)h^2*dh]/[0.6*(π/4)*d^2*(gh)^0.5]=-[h^1.5*dh]/[0.6d^2*(g)^0.5]则263.93-120 =X^2.5/[1.5*d^2*(g)^0.5]以d=0.03m ，g=9.8m*s(-2代入上式得 水深：X第二小题容器为倒葫芦形，比较不规则，比较复杂，不仅要考虑水不断从小孔流出，容器中水的高度是不断变化的，水流速度没有一定的公式，所以要用到微积分解决，还要注意表1的倒葫芦形的不断变化，水深的高度变化是不规则的但仍可以用微积分。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
数学实验报告
1.题目：
某容器盛满水后，底端直径为d0的小孔开启（如图1），根据水力学知识，当水面高度为h时，谁从小孔中流出的速度为v=0.6*（g*h）^0.5（其中g为重力加速度，0.6问哦小孔收缩系数）
1）若容器为倒圆锥形（如图1），现测得容器高和上底直径都为1.2m，小孔直径d为3cm，为水从小孔中流完需要多少时间；2min时水面高度是多少。
2）若容器为倒葫芦形（如图2），现测得容器高1.2m，小孔直径d为3cm，由底端（记x=0）向上每隔0.1m测出容器的直径D（m）如表1所示，问水从小孔中流完需要多少时间；2min时水面高度是多少。
1
表
分析：2.由题知，水从小孔中流出，不仅与容器有关，还与水流速度v=0.6*（2*g*h）^0.5有关。
第一小题容器是圆锥形，比较规则，但是由于水不断从小孔流出，容器中水的高度是不断变化的，水流速度没有一定的公式，所以要用到微积分解决。
由（1）知，水面直径等于水深。
水深为h时，流量为0.6(π/4)d^2*(gh)^0.5，
/4*h^2*dh
π*(d0/2)^2*dt=π0.6*(g*h)^(0.5)*
则水深下降dh所需时间：dt=-[(π/4)h^2*dh]/[0.6(π
/4)d^2*(gh)^0.5]=-[h^1.5*dh]/[0.6d^2*(g)^0.5]
水深由1.2m至0定积分得水从小孔流完的时间：T（其中已知d=0.03m，
g=9.8m*s(-2）
对于第二问：设两分钟（120S）后水深为X m ，由
dt=-[(π/4)h^2*dh]/[0.6*(π/4)*d^2*(gh)^0.5]=-[h^1.5*dh]/[0.6d^2*(g)^0.5] 则263.93-120 =X^2.5/[1.5*d^2*(g)^0.5]
以d=0.03m，g=9.8m*s(-2代入上式得水深：X
第二小题容器为倒葫芦形，比较不规则，比较复杂，不仅要考虑水不断从小孔流出，容器中水的高度是不断变化的，水流速度没有一定的公式，所以要用到微积分解决，还要注意表1的倒葫芦形的不断变化，水深的高度变化是不规则的但仍可以用微积分。
由（2）知容器高1.2m，水深为h时，流量为0.6(π/4)d^2*(gh)^0.5，由于不同高度，倒葫芦形半径不同，用欧拉方程和龙格—库塔方法则水深下降dh所需时间：dt=t(k+1)-t(k)=-[(π/4)h^2*dh]/[0.6(π
/4)d^2*(gh)^0.5]=-[h^1.5*dh]/[0.6d^2*(g)^0.5]
然后利用循环for k=1：length(L)，t(k)=((h(k+1)-h(k))*(π/4)*d(k)^2)/(0.6*(π
/4)*d^2*(g(1.2-h(k)))^0.5)，T=sum(t).可以求得水从小孔流完的总时间。
对于第二问：设两分钟（120S）后水深为X m ，由
S=0,利用条件，当120-s<0.0001时s=s+t(k)，x(k)
把d=0.03m，g=9.8m*s(-2)代入上式得水深：X
相关资料：龙格—库塔方法求解
龙格—库塔方法思想：用[v,v]上若干个点的倒数，对他们做线性组合得到n1n?平均斜率，就可能得到更高阶的精度，这就是龙格—库塔方法思想。
为了演算方便本课程设计采用二阶龙格—库塔公式，即求[v,v]上的二n1n?个倒数，将他们加权平均得平均斜率。
??(v)?f(v,Hh) 我们设有形如：H)v?H(00f(h,v)满足HipscHitz其中函数条件，既满足存在常数H使
f(v,H)?f(v,H)?LH?H2112??v?1??)l(0v去倒数作线性组合和按照下式在nn
??k)(k??HH?h?2n21n?11?k?f(h,v)?1nn??????,1?),y?(?kfvh,?hk0?21nn
（3）
???为待定系数，确定他们的准则是使（3）式有尽量高的精度。
注意到其中,,21
的假设，并对作二元泰勒展开，且利用)vH?H(k nn2???，（3）式可写为：
f,Hf*?Hk?H??f Hv1??k?();k?H(v)?hH2111n2n??;?H(v))?f(v,Hk n1n2?????)h?O(,H(v,H(v))?))hkff(v?)l,H(v?( hk)?Hv(v)?(hfvk?nh1nvnn1nnn22????)hO((H?v(v)?)hH?nn3???????于
是)(h)?v)hH)H?H(v?(??O)hH((v n1nn?12n2故得截断误差为
123???????)(h)H?O?)hH((v)?(?vh?)?TH(v?H(1?)n2n?n?111n?n1221?????,1??就可以使（3）式具有2 阶精度。
容易看出只要：2122
???时的龙格—库塔方以上分析了龙格—库塔方法思想，下面用1?1/2,??21上解决本课程设计题目：法即为改良后的欧拉公式法在MATHAB
模型方程：3..d=0.03m
g=9.8m*s(-2)，h=1.2m，d0=1.2m，开始条件所需dhπ/4)d^2*(gh)^0.5，则
水深下降水深为h时，流量Q为0.6(
/4)d^2*(gh)^0.5]=-[h^1.5*dh]/[0.6d^2*(g)^0.5] πdt=-[(π/4)h^2*dh]/[0.6(时间：），g-9.8m/sT至0定积分得水从小孔流完的时间：（其中已知d=0.03m水深由1.2m X m ，由120S 设两分钟（）后水深为
/4)*d^2*(gh)^0.5]=-[h^1.5*dh]/[0.6d^2*(g)^0.5] /4)h^2*dh]/[0.6*(πdt=-[(π
263.93-120 =X^2.5/[1.5*d^2*(g)^0.5]
则
X 水深：d=0.03，g=9.8代入上式得以
：所需时间dh 用欧拉方程和龙格—库塔方法则水深下降π
/4)h^2*dh]/[0.6(dt=t(k+1)-t(k)=-[(π/4)d^2*(gh)^0.5]=-[h^1.5*dh]/[0.6d^2*(g)^0.5]
k1=0.15*sqrt(g*(x(n)))*d^2/(-43.6359*x(n)^8+213.0457*x(n)^7-414.873*x(n)^6+41 0.2075*x(n)^5-218.8936*x(n)^4+62.553*x(n)^3-8.3215*x(n)^2+0.49619*x(n)+.0148 92)^2;
k2=0.15*sqrt(g*(x(n))-h*k1)*d^2/(-43.6359*(x(n)-h*k1)^8+213.0457*(x(n)-h*k1)^7
-414.873*(x(n)-h*k1)^6+410.2075*(x(n)-h*k1)^5-218.8936*(x(n)-h*k1)^4+62.553*(
x(n)-h*k1)^3-8.3215*(x(n)-h*k1)^2+0.49619*(x(n)-h*k1)+.014892)^2;
x(n+1)=x(n)-h*(k1+k2)/2;
4.编译程序:
(1)
g=9.8;
d=0.03;
syms h
y=-h^(1.5)/(0.6*d^2*sqrt(2*g));
T=int(y,h,1.2,0);
eval(T)
x=((T-120)*(1.5*d^2*sqrt(2*g)))^(0.4);
eval(x)
运行结果如下
>> sy1
ans =
263.9316
ans =
0.9416
(2)
clear;
g=9.8;
d=0.03;
k1=0;
k2=0;
h=0.4;
x(1)=1.2;
for n=1:1000
k1=0.15*sqrt(g*(x(n)))*d^2/(-43.6359*x(n)^8+213.0457*x(n)^7-414.873*x(n)^6+41．0.2075*x(n)^5-218.8936*x(n)^4+62.553*x(n)^3-8.3215*x(n)^2+0.49619*x(n)+.0148 92)^2;
k2=0.15*sqrt(g*(x(n))-h*k1)*d^2/(-43.6359*(x(n)-h*k1)^8+213.0457*(x(n)-h*k1)^7
-414.873*(x(n)-h*k1)^6+410.2075*(x(n)-h*k1)^5-218.8936*(x(n)-h*k1)^4+62.553*(
x(n)-h*k1)^3-8.3215*(x(n)-h*k1)^2+0.49619*(x(n)-h*k1)+.014892)^2;
x(n+1)=x(n)-h*(k1+k2)/2;
end
x(300)
t=0:h:1000*h;
plot(t,x);
axis([0,400,0,1.21]);
运行结果如下
>> sy2
ans =
1.0278
1.2
10.80.6400300350250500100150200点头绪0.94m2min263s水从小孔流完需要，时水面高度是总结：5.。