数学实验课程设计常微分方程数值解
常微分方程数值解实验
有大量的常微分方程虽然从理论上讲,其解是存在的,但我们却无 法求出其解析解,此时,我们需要寻求方程的数值解,在求常微分方程 数值解方面,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的倒葫芦形的不断变化,水深的高度变化是不规则的但仍可以用微积分。
实验报告七常微分方程初值问题的数值解法
浙江大学城市学院实验报告课程名称数值计算方法实验项目名称常微分方程初值问题的数值解法 实验成绩指导老师签名日期2015/12/16 一.实验目的和要求1. 用Matlab 软件掌握求微分方程数值解的欧拉方法和龙格-库塔方法; 2. 通过实例学习用微分方程模型解决简化的实际问题;二.实验内容和原理编程题2-1要求写出Matlab 源程序m 文件,并有适当的注释语句;分析应用题2-2,2-3,2-4,2-5要求将问题的分析过程、Matlab 源程序和运行结果和结果的解释、算法的分析写在实验报告上; 2-1 编程编写用向前欧拉公式和改进欧拉公式求微分方程数值解的Matlab 程序,问题如下:在区间[],a b 内(1)N +个等距点处,逼近下列初值问题的解,并对程序的每一句添上注释语句; Euler 法y=eulera,b,n,y0,f,f1,b1改进Euler 法y=eulerproa,b,n,y0,f,f1,b1 2-2 分析应用题假设等分区间数100n =,用欧拉法和改进欧拉法在区间[0,10]t ∈内求解初值问题()()20(0)10y t y t y '=-⎧⎨=⎩并作出解的曲线图形,同时将方程的解析解也画在同一张图上,并作比较,分析这两种方法的精度; 2-3 分析应用题用以下三种不同的方法求下述微分方程的数值解,取10h = 画出解的图形,与精确值比较并进行分析; 1欧拉法; 2改进欧拉法; 3龙格-库塔方法;2-4 分析应用题考虑一个涉及到社会上与众不同的人的繁衍问题模型;假设在时刻t 单位为年,社会上有人口()x t 人,又假设所有与众不同的人与别的与众不同的人结婚后所生后代也是与众不同的人;而固定比例为r 的所有其他的后代也是与众不同的人;如果对所有人来说出生率假定为常数b ,又如果普通的人和与众不同的人的婚配是任意的,则此问题可以用微分方程表示为:其中变量()()()i p t x t x t =表示在时刻t 社会上与众不同的人的比例,()i x t 表示在时刻t 人口中与众不同的人的数量;1假定(0)0.01,0.02p b ==和0.1r =,当步长为1h =年时,求从0t =到50t =解()p t 的近似值,并作出近似解的曲线图形;2精确求出微分方程的解()p t ,并将你当50t =时在分题b 中得到的结果与此时的精确值进行比较; MATLAB 相关函数求微分方程的解析解及其数值的代入dsolve‘egn1’,‘egn2’,‘x ’ subsexpr,{x,y,…},{x1,y1,…}其中‘egn i ’表示第i 个方程,‘x ’表示微分方程中的自变量,默认时自变量为t ; subs 命令中的expr 、x 、y 为符合型表达式,x 、y 分别用数值x1、x2代入; >>symsxyz>>subs'x+y+z',{x,y,z},{1,2,3} ans= 6>>symsx>>subs'x^2',x,2 ans= 4>>s=dsolve‘12Dy y ∧=+’,‘(0)1y =’,‘x ’ ans= >>symsx >>subss,x,2 ans=右端函数(,)f x y 的自动生成f=inline ‘expr ’,’var1’,‘var2’,……其中’expr ’表示函数的表达式,’var1’,‘var2’表示函数表达式中的变量,运行该函数,生成一个新的函数表达式为fvar1,var2,……; >>f=inline'x+3y','x','y' f=Inlinefunction: fx,y=x+3y >>f2,3 ans= 114,5阶龙格-库塔方法求解微分方程数值解t,x=ode45f,ts,x0,options其中f 是由待解方程写成的m 文件名;x0为函数的初值;t,x 分别为输出的自变量和函数值列向量,t的步长是程序根据误差限自动选定的;若ts=t0,t1,t2,…,tf,则输出在自变量指定值,等步长时用ts=t0:k:tf,输出在等分点;options 用于设定误差限可以缺省,缺省时设定为相对误差310-,绝对误差610-,程序为:options=odeset ‘reltol ’,rt,’abstol ’,at,这里rt,at 分别为设定的相对误差和绝对误差;常用选项见下表;选项名 功能 可选值 省缺值 AbsTol 设定绝对误差正数 RelTol 设定相对误差 正数InitialStep 设定初始步长 正数 自动 MaxStep设定步长上界正数MaxOrder 设定ode15s 的最高阶数 1,2,3,4,5 5 Stats 显示计算成本统计 on,off off BDF 设定ode15s 是否用反向差分on,offoff例:在命令窗口执行>>odefun =inline ‘2*y t y -’,‘t ’,‘y ’;>>[],45(,[0,4],1)t y ode odefun =;ans=>>t y ‘o-’,%解函数图形表示>>45(,[0,4],1)ode odefun %不用输出变量,则直接输出图形 >>[],45(,0:4,1)t y ode odefun =;[],t yans=三.操作方法与实验步骤包括实验数据记录和处理2-1编程编写用向前欧拉公式和改进欧拉公式求微分方程数值解的Matlab 程序,问题如下:在区间[],a b 内(1)N +个等距点处,逼近下列初值问题的解,并对程序的每一句添上注释语句; Euler 法y=eulera,b,n,y0,f,f1,b1改进Euler 法y=eulerproa,b,n,y0,f,f1,b1Euler 法y=eulera,b,n,y0,f,f1,b1 y=zeros1,n+1; y1=y0; h=b-a/n; x=a:h:b; fori=1:n; yi+1=yi+hfxi,yi; end plotx,y holdon%求微分方程的精确解 x1=linspacea,b,100; '精确解为' s=dsolvef1,b1,'x' symsxy1=zeros1,100; for i=1:100y1i=subss,x,x1i; endplotx1,y1,'r'title'红色代表精确解'改进Euler 法y=eulerproa,b,n,y0,f,f1,b1 %求微分方程的数值解 y=zeros1,n+1; y1=y0; h=b-a/n; x=a:h:b; fori=1:n; T1=fxi,yi; T2=fxi+1,yi+hT1; yi+1=yi+h/2T1+T2; end plotx,y holdon%求微分方程的精确解 x1=linspacea,b,100; '精确解为' s=dsolvef1,b1,'x' symsxy1=zeros1,100; fori=1:100 y1i=subss,x,x1i; endplotx1,y1,'r'title'红色代表精确解' 2-2分析应用题假设等分区间数100n =,用欧拉法和改进欧拉法在区间[0,10]t ∈内求解初值问题()()20(0)10y t y t y '=-⎧⎨=⎩并作出解的曲线图形,同时将方程的解析解也画在同一张图上,并作比较,分析这两种方法的精度;1向前欧拉法>>euler0,10,100,10,inline'y-20','x','y','Dy=y-20','y0=10' ans= 精确解为 s= 20-10expx ans= +005Columns1through8(2)改进欧拉法>>eulerpro0,10,100,10,inline'y-20','x','y','Dy=y-20','y0=10' ans= 精确解为 s= 20-10expx ans= +005Columns1through8改进欧拉法的精度比向前欧拉法更高; 2-3分析应用题用以下三种不同的方法求下述微分方程的数值解,取10h = 画出解的图形,与精确值比较并进行分析; 1欧拉法; 2改进欧拉法;2-4分析应用题考虑一个涉及到社会上与众不同的人的繁衍问题模型;假设在时刻t 单位为年,社会上有人口()x t 人,又假设所有与众不同的人与别的与众不同的人结婚后所生后代也是与众不同的人;而固定比例为r 的所有其他的后代也是与众不同的人;如果对所有人来说出生率假定为常数b ,又如果普通的人和与众不同的人的婚配是任意的,则此问题可以用微分方程表示为:其中变量()()()i p t x t x t =表示在时刻t 社会上与众不同的人的比例,()i x t 表示在时刻t 人口中与众不同的人的数量;1假定(0)0.01,0.02p b ==和0.1r =,当步长为1h =年时,求从0t =到50t =解()p t 的近似值,并作出近似解的曲线图形;2精确求出微分方程的解()p t ,并将你当50t =时在分题b 中得到的结果与此时的精确值进行比较;1>>euler0,50,50,,inline'','t','p','Dp=','p0= 1' ans= 精确解为 s=1-99/100expx/500 ans=Columns1through82>>dsolve'Dp=','p0=','t' ans=1-99/100expt/500 >>1-99/100exp ans=与欧拉法求得的精确值差0,0001四.实验结果与分析。
常微分方程数值解法
欧拉方法
总结词
欧拉方法是常微分方程数值解法中最基础的方法之一,其基本思想是通过离散化时间点上的函数值来 逼近微分方程的解。
详细描述
欧拉方法基于微分方程的局部线性化,通过在时间点上逐步逼近微分方程的解,得到一系列离散点上 的近似值。该方法简单易行,但精度较低,适用于求解初值问题。
龙格-库塔方法
总结词
影响
数值解法的稳定性对计算结果的精度和可靠 性有重要影响。
判断方法
通过分析数值解法的迭代公式或离散化方法, 判断其是否具有稳定性和收敛性。
数值解法的收敛性
定义
数值解法的收敛性是指随着迭代次数的增加, 数值解逐渐接近于真实解的性质。
影响
数值解法的收敛性决定了计算结果的精度和 计算效率。
分类
根据收敛速度的快慢,可以分为线性收敛和 超线性收敛等。
判断方法
通过分析数值解法的迭代公式或离散化方法, 判断其是否具有收敛性。
误差分析
定义
误差分析是指对数值解法计算过程中 产生的误差进行定量分析和估计的过 程。
分类
误差可以分为舍入误差、截断误差和 初始误差等。
影响
误差分析对于提高计算精度和改进数 值解法具有重要意义。
分析方法
通过建立误差传递公式或误差估计公 式,对误差进行定量分析和估计。
生物学
生态学、生物种群动态和流行病传播 等问题可以通过常微分方程进行建模
和求解。
化学工程
化学反应动力学、化学工程流程模拟 等领域的问题可以通过常微分方程进 行描述和求解。
经济学
经济系统动态、金融市场模拟和预测 等问题可以通过常微分方程进行建模 和求解。
02 常微分方程的基本概念
常微分方程的定义
常微分方程数值解法课程设计
常微分方程数值解法课程设计常微分方程数值解法课程设计一、背景与意义常微分方程在自然科学、工程技术、社会科学等各个领域都有广泛的应用。
例如,物理学中的牛顿第二定律、电磁学中的麦克斯韦方程、生物学中的种群增长模型等都涉及到常微分方程。
然而,很多常微分方程的解析解很难求得或者不存在,因此数值解法就显得尤为重要。
本次课程设计的目的是使学生掌握常微分方程的数值解法,包括欧拉法、龙格-库塔法等,并能够利用这些方法进行实际问题的建模和计算。
通过本次课程设计,学生将了解数值解法的基本思想、误差分析、稳定性等方面的知识,提高解决实际问题的能力。
二、主要内容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.指导教师评语及成绩:指导教师依据学生的实际报告内容,给出本次实验报告的评价。
常微分方程数值解
淮海工学院实验报告书课程名称:数学实验实验名称:常微分方程数值解班级数学091姓名:耿萍学号:090911107日期:2012.4.6 地点数学实验室指导教师:曹卫平成绩:数理科学系1. 实验目的:1.用MATLAB 软件掌握求微分方程数值解的方法。
2.通过实例学习用微分方程模型解决简化的实际问题。
2. 实验内容:(1)容器盛满水后,底端直径为0d 的小孔开启。
根据水力学知识,当水面高度为h 时,水从小孔中流出的速度为v=0.6 2gh(g 为重力加速度,0.6为空口收缩系数)。
1)若容器为倒圆锥形,现测得容器高和上底面直径均为1.2m ,小孔直径为3cm ,问水从小孔中流完需要多少时间;2分钟时水面高度是多少。
2)若容器为倒葫芦形,现测得容器高1.2m ,小孔直径3cm ,由底端向上每隔0.1m 测出容器的直径如下表所示,问水从小孔中流完需要多少时间;2分钟水面高度是多少。
(图略)(2)一只小船渡过宽为d 的河流,目标是起点A 正对着的另一岸B 点。
已知河水流速为v1与船在静水中的速度v2之比为k 。
1)建立小船航线的方程,求其解析解。
2)设d=100m ,v1=1m/s ,v2=2m/s ,用数值解法求渡河所需x 00.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 11.1 1.2 D0.03 0.05 0.08 0.14 0.19 0.33 0.45 0.68 0.98 1.10 1.20 1.13 1.00时间、任意时刻小船的位置及航行曲线,作图,并与解析解比较。
3.实验步骤:(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.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).可以求得水从小孔流完的总时间。
实验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.了解微分方程的基本概念和求解方法,包括欧拉法、改进的欧拉法和龙格-库塔法。
2. 配置实验环境,准备实验所需的工具和软件,如Python编程语言和相关数值计算库。
3.选择一种微分方程进行研究和求解,可以是一阶、二阶或更高阶的微分方程。
4.使用欧拉法、改进的欧拉法和龙格-库塔法分别求解选定的微分方程,并比较其结果的准确性和稳定性。
5.计算数值解与解析解之间的误差,并进行误差分析和讨论。
6.对比不同数值解法的性能,包括计算时间和计算精度。
7.结果展示和总结,根据实验结果对数值解方法进行评价和选取。
实验结果:以一阶线性常微分方程为例,我们选择经典的“衰减振荡”问题进行实验研究。
该问题的微分方程形式为:dy/dt = -λy其中,λ为正实数。
我们首先使用Python编程语言实现了欧拉法、改进的欧拉法和龙格-库塔法。
进一步,我们选择了λ=0.5和初始条件y(0)=1,使用这三种数值解法求解该微分方程,并比较结果的准确性。
通过对比数值解和解析解可以发现,在短时间内,欧拉法、改进的欧拉法和龙格-库塔法的结果与解析解非常接近。
但随着时间的增加,欧拉法的结果开始偏离解析解,而改进的欧拉法和龙格-库塔法仍然能够提供准确的近似解。
这是因为欧拉法采用线性逼近的方式,误差随着步长的增加而累积,而改进的欧拉法和龙格-库塔法采用更高阶的逼近方式,可以减小误差。
为了更直观地比较不同方法的性能,我们还计算了它们的计算时间。
实验结果显示,欧拉法计算时间最短,而龙格-库塔法计算时间最长。
这表明在计算时间要求较高的情况下,可以选择欧拉法作为数值解方法。
常微分方程的数值解法实验报告
常微分方程的数值解法实验报告实验报告:常微分方程的数值解法摘要:常微分方程(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. 基本概念回顾•微分方程的基本概念与分类•常微分方程的初值问题与边值问题•常微分方程数值解法的分类与相关算法2. 常微分方程数值解法•常微分方程数值解的初步概念与实现方法•常规微分方程数值解法的介绍与实现•常规微分方程数值解法的误差分析3. 偏微分方程数值解法•偏微分方程数值解的基本概念与实现方法•偏微分方程数值解的差分格式与相关算法•偏微分方程数值解的误差分析4. 实例分析•常微分方程实例分析•偏微分方程实例分析•结合实际问题的数值模拟课程设计要求本次课程设计要求学生以实践为主,结合已有数值软件应用数值方法求解微分方程,以期掌握微分方程数值解法的实际应用。
设计过程需包括以下几个步骤:1.客观评价现有数值软件的求解效果。
选取数值软件自带的样例微分方程,按照已有方法求解并进行误差分析。
2.对比不同的数值方法的性能。
选取不同的数值解法,比较不同方法的求解效果,并进行误差分析。
3.解决特殊问题的数值求解。
选取有现实意义的微分方程,根据实际问题需求进行求解并给出数值模拟结果。
4.书面报告。
每个学生需要用Markdown格式书写完整的课程设计报告,包括目标、方法、过程及相关结果和分析等。
课程设计参考书目1.Burden, R. L., Fres, J. D., & Burden, A. M. (2016).Numerical analysis (10th ed.). Boston: Cengage Learning.2.Chapra, S. C., & Canale, R. P. (2011). Numerical methods forengineers (6th ed.). New York: McGraw-Hill.3.Kincd, D., & Cheney, W. (2002). Numerical analysis:Mathematics of scientific computing (3rd ed.). Washington, D.C.:American Mathematical Society.结语微分方程数值解课程设计是一门重要的应用课程,具有一定的挑战性。
常微分方程数值解法
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)中的积分,便导致各种不同的差分法。
- 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.。