数学实验课程设计 常微分方程数值解
常微分方程数值解法
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;
常微分方程数值解实验
有大量的常微分方程虽然从理论上讲，其解是存在的，但我们却无 法求出其解析解，此时，我们需要寻求方程的数值解，在求常微分方程 数值解方面，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日
数学实验——常微分方程数值解
实验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。
求常微分方程的数值解
求常微分方程的数值解一、背景介绍常微分方程（Ordinary Differential Equation，ODE）是描述自然界中变化的数学模型。
常微分方程的解析解往往难以求得，因此需要寻找数值解来近似地描述其行为。
求解常微分方程的数值方法主要有欧拉法、改进欧拉法、龙格-库塔法等。
二、数值方法1. 欧拉法欧拉法是最简单的求解常微分方程的数值方法之一。
它基于导数的定义，将微分方程转化为差分方程，通过迭代计算得到近似解。
欧拉法的公式如下：$$y_{n+1}=y_n+f(t_n,y_n)\Delta t$$其中，$y_n$表示第$n$个时间步长处的函数值，$f(t_n,y_n)$表示在$(t_n,y_n)$处的导数，$\Delta t$表示时间步长。
欧拉法具有易于实现和理解的优点，但精度较低。
2. 改进欧拉法（Heun方法）改进欧拉法又称Heun方法或两步龙格-库塔方法，是对欧拉法进行了精度上提升后得到的一种方法。
它利用两个斜率来近似函数值，并通过加权平均来计算下一个时间步长处的函数值。
改进欧拉法的公式如下：$$k_1=f(t_n,y_n)$$$$k_2=f(t_n+\Delta t,y_n+k_1\Delta t)$$$$y_{n+1}=y_n+\frac{1}{2}(k_1+k_2)\Delta t$$改进欧拉法比欧拉法精度更高，但计算量也更大。
3. 龙格-库塔法（RK4方法）龙格-库塔法是求解常微分方程中最常用的数值方法之一。
它通过计算多个斜率来近似函数值，并通过加权平均来计算下一个时间步长处的函数值。
RK4方法是龙格-库塔法中最常用的一种方法，其公式如下：$$k_1=f(t_n,y_n)$$$$k_2=f(t_n+\frac{\Delta t}{2},y_n+\frac{k_1\Delta t}{2})$$ $$k_3=f(t_n+\frac{\Delta t}{2},y_n+\frac{k_2\Delta t}{2})$$ $$k_4=f(t_n+\Delta t,y_n+k_3\Delta t)$$$$y_{n+1}=y_n+\frac{1}{6}(k_1+2k_2+2k_3+k_4)\Delta t$$三、数值求解步骤对于给定的常微分方程，可以通过以下步骤求解其数值解：1. 确定初值条件：确定$t=0$时刻的函数值$y(0)$。
实验五 常微分方程数值解
F ( t , y ', y '', ⋯ , y ( n ) ) = 0
4
如果未知函数是多元函数，称为偏微分方程。 如果未知函数是多元函数，称为偏微分方程。联系一些未知函数的 一组微分方程称为微分方程组。 一组微分方程称为微分方程组。微分方程中出现的未知函数的导数 的最高阶数称为微分方程的阶。 的最高阶数称为微分方程的阶。若方程中未知函数及其各阶导数都 是一次的，称为线性常微分方程， 是一次的，称为线性常微分方程，一般形式为
y '( t ) = f ( t , y ( t )), t 0 < t < t f y (t0 ) = y 0
y = ( y1 , y2 ,⋯ ym )T , f = ( f1 , f 2 ,⋯ f m )T , y0 = ( y10 , y20 ,⋯ ym 0 )T 其中
所谓数值解法就是寻求y(t ) 在一系列离散节点 t0 < t1 < ⋯ < tn ≤ t f 步长， 上的近似值 y k , k = 0 ,1, ⋯ , n ，称为 hk = t k +1 − t k 步长，通常取为 最简单的数值解法是Euler法。 常量 h 。最简单的数值解法是 法
θ
mlθ '' = mg sin θ ,θ (0) = θ0 ,θ '(0) = 0
问该微分方程是线性的还是非线性的？ 问该微分方程是线性的还是非线性的？是否存 在解析解？如果不存在解析解， 在解析解？如果不存在解析解，能否求出其近 似解？ 似解？
3
【实验准备】 实验准备】
1.微分方程的概念 微分方程的概念 未知的函数以及它的某些阶的导数连同自变量都由一已知 方程联系一起的方程称为微分方程。 方程联系一起的方程称为微分方程。如果未知函数是一元函数 称为常微分方程。 ，称为常微分方程。常微分方程的一般形式为
大学数学实验四_常微分方程数值解
2 小型火箭初始质量为 1400kg，其中包括 1080kg 燃料。火箭竖直向上发射时，燃料燃烧率 为 18kg/s，由此产生 32000N 推力，火箭引擎在燃料用尽时关闭。设火箭上升时空气阻力正 比于速度的平方，比例系数为 0.4kg/m，求引擎关闭瞬间火箭的高度、速度、加速度，及火 箭到达最高点时的高度和加速度，并画出高度、速度、加速度随时间变化的图形。
6
输出结果为
故，当 t =60s 时，火箭的速度、高度、加速度为
v (m/s)
向前欧拉公式
266.94
经典龙格-库塔公式
266.94
ode45(@rocket,t,v)
266.93
h (m) 12131 12131 12131
从上表可以看出，用三种不同方法求得的结果相差不大。
60s 后，火箭作竖直上抛运动，加速度为
4
输出的结果为
输出的火箭的 h – t 图线如下：
经典龙格-库塔公式： n=100000; h=60/n; t=0:h:60; v(1)=0; z(1)=0; g(1)=90/7; for i=1:n
a = v(i); s = z(i); x1 = h*i; k1 = (32000-0.4*a.^2)/(1400-18*x1)-10; x1 = x1 + h/2; b = a + h/2*k1; k2 = (32000-0.4*b.^2)/(1400-18*x1)-10; c = a + h/2*k2; k3 = (32000-0.4*c.^2)/(1400-18*x1)-10; x1 = x1 + h/2;
n = 500 267.0539 n = 100000 266.9432
数学实验课程设计常微分方程数值解
数学实验报告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. 欧拉法 利用matlab中的一个循环语句即可实现欧拉法中的利用第一个节点计
算随后所有节点的工作. function Euler %输出节点的x值和y值, clc %------输入(x0,xn)为求解区间,y0=y(x0)为初始条件,n为区间的 等分个数------x0=0; y0=1; xn=0.5; n=50; %-------------------------------------------------------------------y(1)=y0; h=(xn-x0)/n; x=x0:h:xn; disp('欧拉法结果如下：') disp(['y','(',num2str(x0),')','=',num2str(y0)]); for i=2:n+1 %利用循环语句实现欧拉法利用第一个节 点的函数值计算 y(i)=y(i-1)+h*f(x0,y(i-1)); %随后所有节点函数值的过程. x0=x0+h; disp(['y','(',num2str(x0),')','=',num2str(y(i))]); end plot(x,y,'ro-'); %利用求出的x坐标和y坐标画出解的 近似图形. xlabel('x') ylabel('y') title('欧拉法求出的折线') hold on %-------若原微分方程的理论解能求出，则可画出积分曲线比 较--------------x=0:0.01:0.5; y=x+exp(-x); plot(x,y,'b-') hold off %-----------------------------------------------------------------function y=f(x,z) y=-z+x+1;
常微分方程数值解法课程设计
常微分方程数值解法课程设计常微分方程数值解法课程设计一、背景与意义常微分方程在自然科学、工程技术、社会科学等各个领域都有广泛的应用。
例如，物理学中的牛顿第二定律、电磁学中的麦克斯韦方程、生物学中的种群增长模型等都涉及到常微分方程。
然而，很多常微分方程的解析解很难求得或者不存在，因此数值解法就显得尤为重要。
本次课程设计的目的是使学生掌握常微分方程的数值解法，包括欧拉法、龙格-库塔法等，并能够利用这些方法进行实际问题的建模和计算。
通过本次课程设计，学生将了解数值解法的基本思想、误差分析、稳定性等方面的知识，提高解决实际问题的能力。
二、主要内容1.常微分方程的基本概念：介绍常微分方程的定义、分类、解的存在性和唯一性等基础知识。
2.数值解法的基本思想：介绍数值解法的基本思想，包括离散化、逼近、迭代等，以及数值解法的误差来源和误差估计。
3.欧拉法：介绍欧拉法的基本思想、计算公式、误差分析和稳定性等方面的知识，并通过实例演示欧拉法的应用。
4.龙格-库塔法：介绍龙格-库塔法的基本思想、计算公式、误差分析和稳定性等方面的知识，并通过实例演示龙格-库塔法的应用。
5.实际问题建模与计算：选取实际问题，如物理学中的弹簧振子问题、生物学中的种群增长问题等，利用常微分方程的数值解法进行建模和计算，并对结果进行分析和解释。
三、实施步骤1.理论学习：通过课堂讲解、阅读教材等方式，使学生掌握常微分方程的基本概念、数值解法的基本思想和常用方法。
2.上机实践：安排学生在计算机上利用编程语言实现欧拉法、龙格-库塔法等数值解法，并对简单的常微分方程进行数值计算。
3.实际问题建模与计算：选取实际问题，指导学生利用常微分方程的数值解法进行建模和计算，并对结果进行分析和解释。
4.课程设计报告：要求学生撰写课程设计报告，内容包括问题描述、数学模型、数值解法、计算结果与分析等，以培养学生综合运用所学知识解决实际问题的能力。
四、预期成果通过本次课程设计，学生将能够：1.掌握常微分方程的基本概念和数值解法的基本思想；2.熟练使用欧拉法、龙格-库塔法等常用数值解法；3.能够利用常微分方程的数值解法进行实际问题的建模和计算；4.撰写规范的课程设计报告，具备综合运用所学知识解决实际问题的能力。
数值分析实验报告之常微分方程数值解
数学与计算科学学院实验报告
实验项目名称常微分方程数值解
所属课程名称数值方法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．指导教师评语及成绩：指导教师依据学生的实际报告内容，给出本次实验报告的评价。
实验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); }}。
数值分析教案_常微分方程初值问题数值解法
第九章常微分方程初值问题数值解法图9-1n 作为()n x y 的近似值，得 ()n n y x hf ,)y x ,两边从n x 到1+n x 积分，得()dx x y x f x y x n nx x n n ⎰+=-+1))(,()1 矩形公式计算上式右侧积分，即()()x x x x x d x y x f dx x y x f n nn n⎰⎰++≈11,))(,()n ，得()n n n n y x hf y y ,1+=+，故欧拉法也称为矩形法。
为了达到较高精度的计算公式，对欧拉法进行改进，用梯形公式计算()()([1,2))(,(1++≈+n n n x f x y x f hdx x y x f n 的近似值，得9.2 龙格—库塔法前面讨论的欧拉法与改进的欧拉法都是一步法，即计算y 1+n 时，只用到前一步值。
龙格—库塔(Runge-Kutta)法(简称为R-K 方法)不是通过求导数的方法构造近似公式，而是通过计算不同点上的函数值，并对这些函数值作线性组合，构造近似公式，再把近似公式与解的泰勒展开式进行比较，使前面的若干项相同，从而使近似公式达到一定的阶数。
我们先分析欧拉法与预估—校正法。
对于欧拉法⎩⎨⎧=+=+),(111n n n n y x hf k k y y 每步计算f 的值一次，其截断误差为O (2h )。
对于预估—校正法()()⎪⎪⎩⎪⎪⎨⎧++==++=+121211,,2121k y h x hf k y x hf k k k y y n n n n n n 每步计算f 的值两次，其截断误差为O (3h ).下面对预估—校正法进行改进，将该公式写成更一般的形式()()bh y ah x hf k y x hf k k R k R y y n n n n n n ++==++=+,,2122111 (2.1)其中b a R R ,,,21为待定常数。
选择这些常数的原则是在)(n n x y y =的前提下，使11)(++-n n y x 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、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 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时水面高度是多少。
图1X/m0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2D/ m 0.030.050.080.140.190.330.450.680.981.11.21.131.0表12.分析：由题知，水从小孔中流出，不仅与容器有关，还与水流速度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的倒葫芦形的不断变化，水深的高度变化是不规则的但仍可以用微积分。
由（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 n ,v 1+n ]上若干个点的倒数，对他们做线性组合得到平均斜率，就可能得到更高阶的精度，这就是龙格—库塔方法思想。
为了演算方便本课程设计采用二阶龙格—库塔公式，即求[v n ,v 1+n ]上的二个倒数，将他们加权平均得平均斜率。
我们设有形如： {0)(),()(H v H h v f v H =='其中函数),(v h f 满足HipscHitz 条件，既满足存在常数H 使2121),(),(H H L H v f H v f -≤-按照下式在n v 和)10(≤<+ααl v n 去倒数作线性组合⎪⎩⎪⎨⎧≤<++==++=+1,0),,(),()(12122111βαβαλλhk y h v f k v h f k k k h H H n n nn n n（3）其中αλλ,,21为待定系数，确定他们的准则是使（3）式有尽量高的精度。
注意到)(n n v H H =的假设，并对2k 作二元泰勒展开，且利用H v f f H H f H k *,1+=''='=，（3）式可写为：)()()()())(,())(,()())(,(;))(,();()(22112122111h O v H h v H h O v H v f hk v H v hf v H hk v H l v f k H v H v f k k k h v H H n n n n h n n v n n n n n n n +''+'=+++'=++='==++=+αααααλλ于是)()()()()(32211h O v H h v H h v H H n n n n +''+'++=+αλλλ故得截断误差为)()()21()()1()(32221111h O v H h v H h H v H T n n n n n +''-+'--=-=+++αλλλ容易看出只要： 就可以使（3）式具有2阶精度。
以上分析了龙格—库塔方法思想，下面用1,2/121===αλλ时的龙格—库塔方法即为改良后的欧拉公式法在MATHAB 上解决本课程设计题目：3..模型方程：开始条件h=1.2m ，d0=1.2m ，g=9.8m*s(-2)，d=0.03m水深为h 时，流量Q 为0.6(π/4)d^2*(gh)^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 ） 设两分钟（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.03，g=9.8代入上式得 水深：X用欧拉方程和龙格—库塔方法则水深下降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]k1=0.15*sqrt(g*(x(n)))*d^2/(-43.6359*x(n)^8+213.0457*x(n)^7-414.873*x(n)^6+410.2075*x(n)^5-218.8936*x(n)^4+62.553*x(n)^3-8.3215*x(n)^2+0.49619*x(n)+.014892)^2;1,21212=+=λλαλ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 hy=-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)运行结果如下>> sy1ans =263.9316ans =0.9416(2)clear;g=9.8;d=0.03;k1=0;k2=0;h=0.4;x(1)=1.2;for n=1:1000k1=0.15*sqrt(g*(x(n)))*d^2/(-43.6359*x(n)^8+213.0457*x(n)^7-414.873*x(n)^6+410.2075*x(n)^5-218.8936*x(n)^4+62.553*x(n)^3-8.3215*x(n)^2+0.49619*x(n)+.014892)^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; endx(300)t=0:h:1000*h; plot(t,x);axis([0,400,0,1.21]);运行结果如下 >> sy2ans =1.02780501001502002503003504000.20.40.60.811.2水从小孔流完需要263s ，2min 时水面高度是0.94m 点头绪5.总结：数学实验设计和数学建模有点相似，把实际问题用理论解决，这个题目开始有点兴趣，不过对于第二问时有点难住我，自己把相关知识好好理一下，才有点头绪，在运用的同时更好的掌握知识。