数值分析实验报告之常微分方程数值解
常微分方程数值解实验报告
常微分方程数值解实验报告学院:数学与信息科学专业:信息与计算科学姓名:郑思义学号: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进行计算,并与真解作比较。
数值分析第九章常微分方程数值解法
松弛法
通过迭代更新函数值并逐步放松约束 条件来逼近解,适用于刚性和非刚性 问题。
利用线性组合迭代函数值来逼近解, 具有更高的收敛速度和稳定性。
03
数值解法的稳定性分析
数值解法的稳定性定义
数值解法的稳定性是指当微分方程的初值有微小的扰动时, 其数值解的近似值的变化情况。如果数值解在微小扰动下变 化较小,则称该数值方法是稳定的。
更高的精度和稳定性。
数值逼近法
泰勒级数法
将微分方程的解展开为泰勒级数,通过截断级数来逼 近解。
多项式逼近法
利用多项式来逼近微分方程的解,通过选取合适的基 函数和系数来提高逼近精度。
样条插值法
利用样条函数来逼近微分方程的解,具有更好的光滑 性和连续性。
迭代法
雅可比迭代法
通过迭代更新函数值来逼近微分方程 的解,具有简单易行的优点。
初值和边界条件的处理
根据实际问题,合理设定初值和边界 条件,以获得更准确的数值解。
收敛性和误差分析
对数值解进行收敛性和误差分析,评 估解的精度和稳定性。
数值解法的应用案例分析
人口增长模型
通过数值解法求解人口增长模型,预测未来人口数量,为政策制 定提供依据。
化学反应动力学
利用数值解法研究化学反应的动力学过程,模拟反应过程和结果。
数值分析第九章常微分方 程数值解法
• 引言 • 常微分方程数值解法的基本思想 • 数值解法的稳定性分析 • 数值解法的收敛性和误差分析 • 数值解法的实现和应用案例
01
引言
常微分方程的应用背景
自然科学
描述物理、化学、生物等自然 现象的变化规律。
工程领域
控制系统设计、航天器轨道计 算等。
实验报告——常微分方程的数值解法
实验报告实验项目名称常微分方程的数值解法实验室数学实验室所属课程名称微分方程数值解实验类型上机实验实验日期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)(),(y a y bx a y x f dxdy, (9-1) 的数值解法具有典型性。
常微分方程的解能用初等函数、特殊函数或它们的级数与积分表达的很少。
用解析方法只能求出线性常系数等特殊类型的方程的解。
对非线性方程来说,解析方法一般是无能为力的,即使某些解具有解析表达式,这个表达式也可能非常复杂而不便计算。
因此研究微分方程的数值解法是非常必要的。
只有保证问题(9-1)的解存在唯一的前提下,研究其数值解法或者说寻求其数值解才有意义。
由常微分方程的理论知,如果(9-1)中的),(y x f 满足条件(1)),(y x f 在区域} ),({+∞<<∞-≤≤=y b x a y x D ,上连续; (2)),(y x f 在上关于满足Lipschitz 条件,即存在常数,使得y y L y x f y x f -≤-),(),(则初值问题(9-1)在区间],[b a 上存在惟一的连续解)(x y y =。
在下面的讨论中,我们总假定方程满足以上两个条件。
所谓数值解法,就是求问题(9-1)的解)(x y y =在若干点b x x x x a N =<<<<= 210处的近似值),,2,1(N n y n =的方法。
),,2,1(N n y n =称为问题(9-1)的数值解,n n x x h -=+1称为由到1+n x 的步长。
今后如无特别说明,我们总假定步长为常量。
建立数值解法,首先要将微分方程离散化,一般采用以下几种方法: (1) 用差商近似导数在问题(9-1)中,若用向前差商hx y x y n n )()(1-+代替)(n x y ',则得)1,,1,0( ))(,()()(1-=≈-+N n x y x f hx y x y n n n n n)(n x y 用其近似值代替,所得结果作为)(1+n x y 的近似值,记为1+n y ,则有 1(,) (0,1,,1)n n n n y y hf x y n N +=+=-这样,问题(9-1)的近似解可通过求解下述问题100(,) (0,1,,1)()n n n n y y hf x y n N y y x +=+=-⎧⎨=⎩(9-2)得到,按式(9-2)由初值经过步迭代,可逐次算出N y y y ,,21。
数值分析常微分方程数值解法
第8页/共105页
➢ 数值积分方法(Euler公式)
设将方程 y=f (x, y)的两端从 xn 到xn+1 求积分, 得
y( xn1) y( xn )
xn1 f ( x, y( x))dx :
xn
xn1 F ( x)dx
xn
用不同的数值积分方法近似上式右端积分, 可以得到计算 y(xn+1)的不同的差分格 式.
h2 2
y''( )
Rn1
:
y( xn1)
yn1
h2 2
y''( )
h2 2
y''( xn ) O(h3 ).
局部截断误差主项
19
第20页/共105页
➢ 向后Euler法的局部截断误差
向后Euler法的计算公式
yn1 yn hf ( xn1, yn1 ), n 0, 1, 2,
定义其局部截断误差为
y 计算 的n递1 推公式,此类计算格式统称为差分格式.
3
第4页/共105页
数值求解一阶常微分方程初值问题
y' f ( x, y), a x b,
y(a)
y0
难点: 如何离散 y ?
➢ 常见离散方法
差商近似导数 数值积分方法 Taylor展开方法
4
第5页/共105页
➢ 差商近似导数(Euler公式)
(0 x 1)
y(0) 1.
解 计算公式为
yn1
yn
hfn
yn
h( yn
2xn ), yn
y0 1.0
n 0, 1, 2,
取步长h=0.1, 计算结果见下表
13
常微分方程数值解实验报告
常微分方程数值解实验报告实验报告:常微分方程数值解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日
求常微分方程的数值解
求常微分方程的数值解一、背景介绍常微分方程(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)$。
大学数学实验四_常微分方程数值解
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
微分方程数值解法实验报告
微分方程数值解法实验报告2篇微分方程数值解法实验报告(一)在实际科学与工程问题中,我们经常会遇到微分方程的求解。
然而,许多微分方程往往没有解析解,这就需要我们利用数值方法来获得近似解。
本篇实验报告将介绍两种常见的微分方程数值解法:欧拉方法和改进的欧拉方法。
一、欧拉方法欧拉方法是最简单的微分方程数值解法之一。
其基本原理为离散化微分方程,将微分方程中的导数用差商代替。
设要求解的微分方程为dy/dx = f(x, y),步长为h,则可用以下公式进行递推计算:y_{n+1} = y_n +hf(x_n, y_n)二、算法实现为了对欧拉方法进行数值实验,我们以一阶线性常微分方程为例:dy/dx = x - y, y(0) = 1步骤如下:(1)选择合适的步长h和求解区间[a, b],这里我们取h=0.1,[a, b] = [0, 1];(2)初始化y_0 = 1;(3)利用欧拉方法递推计算y_{n+1} = y_n + 0.1(x_n - y_n);(4)重复步骤(3),直到x_n = 1。
三、实验结果与讨论我们通过上述步骤得到了在[0, 1]上的近似解y(x)。
下图展示了欧拉方法求解的结果。
从图中可以看出,欧拉方法得到的近似解与精确解有一定的偏差。
这是因为欧拉方法只是通过递推计算得到的近似解,并没有考虑到更高阶的误差。
所以在需要高精度解时,欧拉方法并不理想。
四、改进的欧拉方法针对欧拉方法的不足,我们可以考虑使用改进的欧拉方法(也称为改进的欧拉-柯西方法)。
它是通过利用前后两个步长欧拉方法得到的结果作为差商的中间项,从而提高了求解精度。
一阶线性常微分方程的改进欧拉方法可以表示为:y_{n+1} = y_n + h(f(x_n, y_n) + f(x_{n+1}, y_n + hf(x_n,y_n)))/2五、算法实现与结果展示将改进的欧拉方法应用于前述的一阶线性常微分方程,我们同样选择h=0.1,[a, b] = [0, 1]。
《数值分析》第5讲:常微分方程数值解
1 2
y(3) ( xn )2 p2h3
o(h4 )
对照
y( xn1 )
yn
y'( xn )h
1 2
y''( xn )h2
第五章:常微分方程数值解
得
212 p2
1
1
可解得
yn1 yn h(1K1
K1 f ( xn , yn )
2
K
2
)
K
2
f ( xn p , yn
phK1 )
12
第五章:常微分方程数值解
§5.3 Lunge-Kutta方法
依据精度要求的待定系数法
x xn
x n p n1
1、二阶Lunge-Kutta方法(P114-P116)
令 xn p xn ph( xn, xn1), 0 p 1 yn1 yn h(1K1 2 K2 ) 加权平均斜率
K1 f ( xn , yn ) y(x) 在点 ( xn , 的yn斜) 率
例 P107
初值问题
y' y 2x / y
y(0)
1
0 x 1
补 充 : 一 阶 线 性 方 程dy P( x) y Q( x)的 解
dx
常数变易法
y e p( x)dx
Q( x)e
p( x )dx
dx
C
非线性方程 dy P( x) y Q( x) yn Bernoulli型方程 dx
yn1 yn h(1K1 2 K 2 )
yn h[1 y'( xn ) 2 y'( xn p )]
yn
h1 y'( xn ) 2h
1 y'( xn ) 1
数值分析实验报告之常微分方程数值解
数学与计算科学学院实验报告
实验项目名称常微分方程数值解
所属课程名称数值方法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方法和其他变体中具有不同的取值。
数值分析第九章常微分方程数值解法
数值分析第九章常微分方程数值解法常微分方程数值解法是数值分析中非常重要的一部分内容。
常微分方程是描述自然现象中动态变化规律的数学模型,解常微分方程可以揭示系统的变化趋势和规律。
然而,大多数常微分方程是无法通过解析方法求出解的,因此需要借助计算机进行数值计算。
数值解常微分方程方法主要包括:Euler方法、改进的Euler方法、四阶Runge-Kutta方法和龙格-库塔方法。
Euler方法是最简单的方法之一,它采用的是一阶Taylor展开式。
将待求的函数值与函数的一阶导数值代入Taylor展开式中,可以得到函数值在下一个时间步长上的近似值。
Euler方法的优点是简单易于实现,但其精度不够高,容易积累误差。
改进的Euler方法是对Euler方法的改进,它通过使用中间点上的导数值来减小误差。
改进的Euler方法的精度相比Euler方法要高一些,但仍然不够高。
四阶Runge-Kutta方法是目前使用较为广泛的数值解常微分方程的方法之一、它通过计算不同时间点上的斜率来估计函数值,在多个时间点上计算斜率的平均值来提高精度。
四阶Runge-Kutta方法的精度比Euler方法和改进的Euler方法要高,但计算量也相对较大。
龙格-库塔方法是数值解常微分方程中最常用的方法之一、它是四阶Runge-Kutta方法的延伸,通过计算不同时间点上的斜率来估计函数值,然后利用这些估计值计算更准确的斜率,在不同步长上进行迭代计算,直到满足所需精度。
龙格-库塔方法的精度比四阶Runge-Kutta方法要高,但计算量也相对较大。
除了以上几种方法外,还有一些其他数值解常微分方程的方法,如Adams法、Gear法等。
这些方法在不同场景下有着不同的适用性和优劣势。
总结起来,数值解常微分方程方法是研究常微分方程数值计算中的重要内容。
不同的方法有着不同的精度和计算量,可以根据具体问题的特点选择合适的方法进行数值计算。
然而,需要注意的是,数值解只是在给定的步长下对函数的近似值,可能会引入误差。
数学建模-- 常微分方程数值解及实验
注意:
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 ))]
数值分析21(常微分方程数值解)
下面以2级方法为例子具体介绍龙格-库塔法。 yn1 yn hc1 K 1 hc2 K 2
K 1 f ( x n , yn ) K 2 f ( xn 2 h, yn h21 K 1 )
构造的基本思想是选择适当的系数使得方法的局部 截断误差阶数尽可能高。
Tn1 y( xn1 ) yn hc1 f ( xn , yn ) hc2 f ( xn 2h, yn h21 f ( xn , yn ))
y ( xn1 ) y ( xn ) h
f ( , y( )), 其中 [ xn , xn1 ]
h yn1 yn [ f ( xn , yn ) f ( xn1 , yn1 )] 2 是否可以推广改进的Euler方法?
(Runge-Kutta)龙格-库塔法
RK方法是一大类的方法, 其基本思想是采用如下形式:
yn1 yn hf ( xn1 , yn1 ( xn1 )) 隐式Euler法
梯形公式:
x n 1 xn
h f ( x , y( x ))dx [ f ( xn , yn ) f ( xn1 , y( xn1 ))] 2
20:24
h yn1 yn [ f ( xn , yn ) f ( xn1 , y( xn1 ))] 2
c1 c2 1, c22 1 / 2, c221 1 / 2
故方程组有无穷多组解, 每一组来自构成的RK公式的阶数 都是2, 都叫做二阶RK公式。
c1 c2 1 / 2, 2 21 1就是改进Euler法。
20:24
16/27
三阶Range-Kutta公式一般形式
《数值分析》 23
数值分析2.5__常微分方程的数值解
三、单步法的局部截断误差和精度 有关) 单步法的一般形式为 : (ϕ与f 有关) 显式单步法 + hϕ ( x n , y n , y n +1 , h)
yn +1 = yn + hϕ ( xn , yn , h)
整体截断误差: 整体截断误差:从x0开始,考虑每一步产生的 开始, 误差,直到x 误差,直到 n,则有误差 e = y ( x ) − y
y′′ = f ( x, y, y′) a ≤ x ≤ b (3) y ( a ) = y0 , y (b) = yn
(1),(2)式称为初值问题,(3)式称为边值问题。 , 式称为初值问题 式称为初值问题, 式称为边值问题 式称为边值问题。 在实际应用中还经常需要求解常微分方程组: 在实际应用中还经常需要求解常微分方程组:
xn 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 y(xn) 0 0.1923 0.3448 0.4412 0.4878 0.5000 0.4918 0.4730 yn 0 0.2000 0.3840 0.5170 0.5824 0.5924 0.5705 0.5354 yn-y(xn) 0 0.0077 0.0392 0.0758 0.0946 0.0924 0.0787 0.0624
并将数值解和该问题的解析解比较。 并将数值解和该问题的解析解比较。 x 解析解: y ( x) = 2 1+ x 方法的具体格式: 解:Euler方法的具体格式: 方法的具体格式
y n +1
yn 2 = y n + h( − 2 y n ) xn
取h=0.2, xn=nh,(n=0,1,2…,15), f(x,y)=y/x – 2y2 计算中取f(0,0)=1. 计算结果如下: 计算结果如下: 计算中取
实验八 常微分方程初值问题数值解法报告
实验八 常微分方程初值问题数值解法一、基本题科学计算中经常遇到微分方程(组)初值问题,需要利用Euler 法,改进Euler 法,Rung-Kutta 方法求其数值解,诸如以下问题:(1) ()⎪⎩⎪⎨⎧=-='004y xy y x y 20≤<x分别取h=0.1,0.2,0.4时数值解。
初值问题的精确解245x y e -=+。
(2) ()⎩⎨⎧=--='0122y y x y 01≤≤-x用r=3的Adams 显式和预 - 校式求解取步长h=0.1,用四阶标准R-K 方法求值。
(3)()()()100010321331221==-='⎪⎩⎪⎨⎧-='-='='y y y y y y y y y 10≤≤x用改进Euler 法或四阶标准R-K 方法求解取步长0.01,计算(0.05),(0.1y y y 数值解,参考结果 123(0.15)0.9880787,(0.15)0.1493359,(0.15)0.8613125y y y ≈-≈≈。
(4)利用四阶标准R- K 方法求二阶方程初值问题的数值解(I )()()⎩⎨⎧='==+'-''10,00023y y y y y 02.0,10=≤≤h x(II)()()()⎩⎨⎧='==+'--''00,10011.02y y y y y y 1.0,10=≤≤h x(III)()()⎪⎩⎪⎨⎧='=+='00,101y y e y y x 1.0,20=≤≤h x(IV)()()⎩⎨⎧='==+''00,100sin y y y y 2.0,40=≤≤h x二、应用题1. 小型火箭初始质量为900千克,其中包括600千克燃料。
火箭竖直向上发射时燃料以15千克/秒的速率燃烧掉,由此产生30000牛顿的恒定推力.当燃料用尽时引擎关闭。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
数学与计算科学学院实验报告
实验项目名称常微分方程数值解
所属课程名称数值方法B
实验类型验证
实验日期 2013.11.11
班级
学号
姓名
成绩
【实验过程】(实验步骤、记录、数据、分析)
注:以下图形是通过Excel 表格处理数据得出,并未通过MATLAB 编程序所得!
1、1(0)1dy
y x dx y ⎧=-++⎪⎨⎪=⎩
由题可知精确解为:x y x e -=+,当x=0时,y(x)=0。
h=0.1
表1 h=0.1时三个方法与精确值的真值表
图1 h=0.1时三个方法走势图
步长 Euler 法 预估校正法 经典四阶库 精确值 0.1 1.010000 1.005000 1.004838 1.249080 0.2 1.029000 1.019025 1.018731 1.055455 0.3 1.056100 1.041218 1.040818 1.091217 0.4 1.090490 1.070802 1.070320 1.131803 0.5 1.131441 1.107076 1.106531 1.176851 0.6 1.178297 1.149404 1.148812 1.226025 0.7 1.230467 1.197211 1.196586 1.279016 0.8 1.287421 1.249975 1.249329 1.335536 0.9 1.348678 1.307228 1.306570 1.395322 1.0
1.413811
1.368541
1.367880
1.458127
h=0.05(此时将源程序中i的围进行扩大,即for(i=0;i<20;i++))
表2 h=0.05时三个方法与精确值的真值表步长Euler法预估校正法经典四阶库精确值
0.05 1.002500 1.001250 1.001229 1.011721
0.10 1.007375 1.004877 1.004837 1.024908
0.15 1.014506 1.010764 1.010708 1.039504
0.20 1.023781 1.018802 1.018731 1.055455
0.25 1.035092 1.028885 1.028801 1.072710
0.30 1.048337 1.040915 1.040818 1.091217
0.35 1.063421 1.054795 1.054688 1.110931
0.40 1.080250 1.070436 1.070320 1.131801
0.45 1.098737 1.087752 1.087628 1.153791
0.50 1.118800 1.106662 1.106531 1.176851
0.55 1.140360 1.127087 1.126950 1.200942
0.60 1.163342 1.148954 1.148812 1.226025
0.65 1.187675 1.172193 1.172046 1.252062
0.70 1.213291 1.196736 1.196585 1.279016
0.75 1.240127 1.222520 1.222367 1.306852
0.80 1.268121 1.249485 1.249329 1.335536
0.85 1.297215 1.277572 1.277415 1.365037
0.90 1.327354 1.306728 1.306570 1.395322
0.95 1.358486 1.336900 1.336741 1.426362
1.00 1.390562 1.368039 1.367880 1.458127
图2 h=0.05时三个方法走势图
2、(0)1x
dy xe y dx y -⎧=-⎪⎨⎪=⎩
由题可知精确解为:21
(2)2
x x e -+,当x=0时,y(x)=0。
h=0.1
表3 h=0.1时三个方法与精确值的真值表
步长 Euler 法 预估校正法 经典四阶库 精确值 0.1 0.900000 0.909625 0.909428 0.929533 0.2 0.819249 0.835927 0.835593 0.872564 0.3 0.754433 0.776081 0.775655 0.826822 0.4 0.702726 0.727671 0.727189 0.790348 0.5 0.661726 0.688636 0.688127 0.761457 0.6 0.629396 0.657225 0.656711 0.738709 0.7 0.604018 0.631957 0.631453 0.720874 0.8 0.584147 0.611582 0.611100 0.706908 0.9 0.568575 0.595050 0.594599 0.695927 1.0 0.556297 0.581487 0.581072
0.687191
图3 h=0.1时三个方法走势图
h=0.05(此时将源程序中i的围进行扩大,即for(i=0;i<20;i++))
表4 h=0.05时三个方法与精确值的真值表
步长Euler法预估校正法经典四阶库精确值0.050.9500000.9524520.9524270.962924 0.100.9049040.9094740.9094280.929533 0.150.8642840.8706700.8706060.899511 0.200.8277410.8356710.8355920.872564 0.250.7949080.8041370.8040470.848419 0.300.7654470.7757550.7756550.826822 0.350.7390430.7502320.7501250.807538 0.40 0.7154070.7273020.7271890.790348 0.45 0.6942720.7067150.7065990.775050 0.50 0.6753940.6882450.6881260.761457 0.55 0.6585460.6716820.6715610.749397 0.60 0.6435190.6568300.6567100.738709 0.65 0.6301240.6435140.6433950.729247 0.70 0.6181850.6315700.6314530.720874 0.75 0.6075410.6208480.6207330.713466 0.80 0.5980460.6112110.6111000.706908 0.85 0.5895650.6025350.6024260.701094 0.90 0.5819760.5947030.5945990.695927
0.95 0.5751670.5876120.5875120.691320
1.00 0.5690350.5811670.5810710.687191
图4 h=0.05时三个方法走势图
附录1:源程序
附录2:实验报告填写说明
1.实验项目名称:要求与实验教学大纲一致。
2.实验目的:目的要明确,要抓住重点,符合实验教学大纲要求。
3.实验原理:简要说明本实验项目所涉及的理论知识。
4.实验环境:实验用的软、硬件环境。
5.实验方案(思路、步骤和方法等):这是实验报告极其重要的容。
概括整个实验过程。
对于验证性实验,要写明依据何种原理、操作方法进行实验,要写明需要经过哪几个步骤来实现其操作。
对于设计性和综合性实验,在上述容基础上还应该画出流程图、设计思路和设计方法,再配以相应的文字说明。
对于创新性实验,还应注明其创新点、特色。
6.实验过程(实验中涉及的记录、数据、分析):写明具体实验方案的具体实施步骤,包括实验过程中的记录、数据和相应的分析。
7.实验结论(结果):根据实验过程中得到的结果,做出结论。
8.实验小结:本次实验心得体会、思考和建议。
9.指导教师评语及成绩:指导教师依据学生的实际报告容,给出本次实验报告的评价。