龙格-库塔求解微分方程
龙格库塔法解微分方程组
龙格库塔法解微分方程组引言微分方程组是数学中经常遇到的问题，在物理、工程和自然科学中都有广泛的应用。
为了求解微分方程组，我们需要利用数值方法来逼近解析解。
本文将介绍一种常用的数值方法——龙格库塔法（Runge-Kutta method），并探讨如何利用该方法来解微分方程组。
龙格库塔法概述龙格库塔法是一种迭代数值方法，用于求解常微分方程的初值问题。
它的主要思想是将微分方程的解进行离散化，将其转化为一系列的逼近值。
龙格库塔法的基本步骤如下：1.确定步长h和迭代次数n。
2.初始化初始条件，并假设第一个逼近值为y(xi)。
3.依次计算每个逼近值，直到得到y(xi+n*h)为止。
在每个迭代步骤中，龙格库塔法根据前一步的逼近值来计算下一步的逼近值。
该方法具有高阶精度和较好的稳定性，在实际应用中广泛使用。
单一微分方程的龙格库塔法首先，我们来看如何利用龙格库塔法来解一阶常微分方程。
以方程dy/dx = f(x, y)为例，其中f(x, y)为给定的函数。
步骤一：确定步长和迭代次数选择合适的步长h和迭代次数n来进行数值计算。
步长h决定了离散化的精度，而迭代次数n决定了逼近解的数目。
步骤二：初始化条件并计算逼近值设初始条件为y(x0) = y0，其中x0为起始点，y0为起始点处的函数值。
我们先通过欧拉法计算出y(x0 + h)的逼近值，然后再通过该逼近值来计算下一个逼近值。
逼近值的计算公式如下：k1 = h * f(x0, y0)k2 = h * f(x0 + h/2, y0 + k1/2)k3 = h * f(x0 + h/2, y0 + k2/2)k4 = h * f(x0 + h, y0 + k3)y(x0 + h) = y0 + 1/6 * (k1 + 2k2 + 2k3 + k4)步骤三：重复步骤二直到得到y(xi+n*h)依次利用上一步计算出的逼近值来计算下一个逼近值，直到得到y(xi+n*h)为止。
微分方程组的龙格库塔法对于一阶微分方程组的初值问题，我们可以将其转化为向量形式。
matlab龙格库塔方法求解二元二阶常微分方程组
matlab龙格库塔方法求解二元二阶常微分方程组文章标题：深入探讨matlab中的龙格库塔方法及其在求解二元二阶常微分方程组中的应用摘要：在科学与工程领域，常常需要求解复杂的微分方程组，而matlab作为一种强大的数学工具，提供了许多求解微分方程组的方法。
本文将深入探讨matlab中的龙格库塔方法及其在求解二元二阶常微分方程组中的应用，以便读者全面理解该方法并能灵活应用于实际问题中。
正文：一、介绍龙格库塔方法龙格-库塔法（Runge-Kutta methods）是一种数值求解常微分方程的方法，通过将微分方程的解进行离散化，将微分方程转化为差分方程，从而进行数值求解。
龙格库塔方法通过迭代计算，能够得到微分方程的数值解，广泛应用于科学计算和工程技术领域。
二、matlab中的龙格库塔方法在matlab中，龙格库塔方法通过ode45函数实现，该函数能够对一阶或高阶常微分方程进行数值求解。
用户可以通过设定初始条件、微分方程表达式，以及积分区间等参数，快速得到微分方程的数值解。
ode45函数采用自适应步长的方式进行求解，能够有效解决微分方程解的数值稳定性和精确度问题。
三、龙格库塔方法在求解二元二阶常微分方程组中的应用考虑如下形式的二元二阶常微分方程组：$$\begin{cases}y_1' = f_1(t, y_1, y_2) \\y_2' = f_2(t, y_1, y_2)\end{cases}$$其中，$y_1(t)$和$y_2(t)$是未知函数，$f_1(t, y_1, y_2)$和$f_2(t,y_1, y_2)$分别表示其对应的函数表达式。
通过matlab中的ode45函数，可以将该二元二阶常微分方程组转化为一阶常微分方程组的形式，然后利用龙格库塔方法进行数值求解。
设定初始条件$y_1(0) = y1_0, y_2(0) = y2_0$，对应的一阶方程组为：$$\begin{cases}u_1' = u_3 \\u_2' = u_4 \\u_3' = f_1(t, u_1, u_2) \\u_4' = f_2(t, u_1, u_2)\end{cases}$$其中，$u_1(t) = y_1(t), u_2(t) = y_2(t), u_3(t) = y_1'(t), u_4(t) =y_2'(t)$，通过ode45函数求解该一阶常微分方程组即可得到原二元二阶常微分方程组的数值解。
龙格-库塔法
四阶龙格－库塔法求解常微分方程的初值问题1.算法原理对于一阶常微分方程组的初值问题⎪⎪⎪⎩⎪⎪⎪⎨⎧=⋯⋯==⋯⋯=⋯⋯⋯⋯=⋯⋯=0020********'212'2211'1)(,,)(,)())(,),(),(,()())(,),(),(,()())(,),(),(,()(n n n n n n n y x y y x y y x y x y x y x y x f x y x y x y x y x f x y x y x y x y x f x y ， 其中b x a ≤≤。
若记Tn Tn Tn y x f y x f y x f y x f y y y y x y x y x y y x y )),(,),,(),,((),(),,,())(),(),(()(2102010021⋯⋯=⋯⋯=⋯⋯=，，则可将微分方程组写成向量形式⎩⎨⎧=≤≤=0')()),(,()(y a y b x a x y x f x y微分方程组初值问题在形式上和单个微分方程处置问题完全相同，只是数量函数在此变成了向量函数。
因此建立的单个一阶微分方程初值问题的数值解法，可以完全平移到求解一阶微分方程组的初值问题中，只不过是将单个方程中的函数转向向量函数即可。
标准4阶R-K 法的向量形式如下：⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧++=++=++==++++=+),()21,2()21,2(),()22(61342312143211K y h x hf K K y h x hf K K y h x hf K y x hf K K K K K y y n n n n n n n n n n 其分量形式为n j K y K y K y h x hf K K y K y K y h x hf K K y K y K y h x hf K y y y x hf K K K K K y y n ni i i i j j n nii i i j j n nii i i j j ni i i i j j j j j j i j i j ,,2,1).,,,;(),2,2,2;2(),2,2,2;2(),,,,;(),22(6132321314222212131212111221143211,1,⋯⋯=⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧+⋯⋯+++=+⋯⋯+++=+⋯⋯+++=⋯⋯=++++=++，，2.程序框图3.源代码%该函数为四阶龙格-库塔法function [x,y]=method(df,xspan,y0,h)%df为常微分方程，xspan为取值区间，y0为初值向量，h为步长x=xspan(1):h:xspan(2);m=length(y0);n=length(x);y=zeros(m,n);y(:,1)=y0(:);for i=1:n-1k1=feval(df,x(i),y(:,i));k2=feval(df,x(i)+h/2,y(:,i)+h*k1/2);k3=feval(df,x(i)+h/2,y(:,i)+h*k2/2);k4=feval(df,x(i)+h,y(:,i)+h*k3);y(:,i+1)=y(:,i)+h*(k1+2*k2+2*k3+k4)/6;end%习题9.2clear;xspan=[0,1];%取值区间h=0.05;%步长y0=[-1,3,2];%初值df=@(x,y)[y(2);y(3);y(3)+y(2)-y(1)+2*x-3];[xt,y]=method(df,xspan,y0,h)syms t;yp=t*exp(t)+2*t-1;%微分方程的解析解yp1=xt.*exp(xt)+2*xt-1%计算区间内取值点上的精确解[xt',y(1,:)',yp1']%y(1,:)为数值解，yp1为精确解ezplot(yp,[0,1]);%画出解析解的图像hold on;plot(xt,y(1,:),'r');%画出数值解的图像4.计算结果。
四阶龙格库塔法(Runge-Kutta)求解微分方程
四阶龙格库塔法（Runge-Kutta ）求解微分方程张晓颖（天津大学 材料学院 学号：1012208027）1 引言计算传热学中通常需要求解常微分方程。
这类问题的简单形式如下：{),(')(00y x f y y x y == （1）虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些特殊类型的方程，实际问题中的多数微分方程需要采用数值解法求解。
初值问题（1）的数值解法有个基本特点，它们采取“步进式”，即求解过程顺着节点排序一步一步向前推进。
这类算法是要给出用已知信息y n 、 y n −i ……计算y n +1的递推公式即可。
2 龙格库塔法（Runge-Kutta ）介绍假设对于初值问题（1）有解 y = y (x ) ，用 Taylor 展开有：......)(!3)(!2)()()(321+'''+''+'+=+n n n n n x y h x y h x y h x y x y （2） 龙格库塔法（Runge-Kutta ）实质上是间接的使用 Taylor 级数法的一种方法。
对于差商hx y x y n n )()(1-+，根据微分中值定理，存在 0 < θ < 1 ，使得：)()()(1h x y hx y x y n n θ+'=-+ （3）于是对于 y = y (x ) ，可得到：))(,()()(1h x y h x hf x y x y n n n n θθ+++=+ （4）设))(,(*h x y h x f K n n θθ++=，为区间 [x n , x n +1 ] 上的平均斜率。
四阶龙格库塔格式中的*K 由下式计算得到：⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎨⎧++=++=++==++++=+),()2,2()2,2(),()22(6342312143211hK 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 n n n n nn n n n n （5） 四阶龙格库塔法（Runge-Kutta ）的每一步需要四次计算函数值f ，其截断误差更低，计算的精度更高。
计算方法 15 龙格库塔-常微分方程
以上公式就是三阶龙格-库塔公式
计算方法（2016/2017 第一学期）
西南科技大学
制造科学与工程学院
11
四阶龙格-库塔公式
四阶龙格-库塔公式：
h yn1 yn 6 ( k1 2k2 2k3 k4 ) k1 f ( xn , yn ) h hk1 ) k2 f ( xn , yn 2 2 h hk2 k3 f ( xn 2 , yn 2 ) k4 f ( xn h, yn hk3 )
以上公式就是四阶龙格-库塔公式
也称作经典龙格-库塔公式
计算方法（2016/2017 第一学期） 西南科技大学 制造科学与工程学院
12
高阶龙格-库塔公式
高阶龙格-库塔公式：
yn1 yn h( λ1k1 λ2 k2 λ3 k3 λm km ) k1 f ( xn , yn ) k2 f ( xn 2 h, yn 21hk1 ) k3 f ( xn 3 h, yn 31hk1 32 hk2 ) km f ( xn m h, yn m1hk1 m ( m 1) hkm 1 ) i (i 2,, m) 和 ij (i 2,, m; 其中 i (i 1,, m)，
14
龙格-库塔公式
由于龙格-库塔法的导出基于泰勒展开，故精度
主要受解函数的光滑性影响。
对于光滑性不好的解，最好采用低阶算法来求解，
而将步长 h 取小。
计算方法（2016/2017 第一学期）
西南科技大学
制造科学与工程学院
15
例题
例：利用四阶龙格-库塔方法求解微分方程的
matlab迭龙格库塔法解常微分方程
一、介绍迭龙格-库塔法（Runge-Kutta method）是一种数值求解常微分方程（ODE）的常用方法。
它是由卡尔·迭龙格（Carl Runge）和马丁·威尔黑尔姆·库塔（Wilhelm Kutta）在20世纪初提出的，该方法以两位数值分析家的名字来命名。
二、简单描述迭龙格-库塔法是通过数值逼近的方式，来计算常微分方程的近似解。
它是一种显式求解方法，适用于解非线性常微分方程和具有较大阶数的常微分方程。
三、数学原理迭龙格-库塔法主要是通过将微分方程转化为差分方程，利用数值解的方式来逼近微分方程的解。
它是一种显式方法，通过不断迭代得到下一个时间步的近似解。
四、matlab中的应用在matlab中，可以使用ode45函数来调用迭龙格-库塔法求解常微分方程。
ode45函数是matlab中集成的一个函数，通过调用ode45函数，可以直接求解常微分方程的数值解。
五、实例演示下面通过一个简单的例子来演示如何使用matlab中的ode45函数来求解常微分方程。
我们考虑一个简单的一阶常微分方程：dy/dt = -y初始条件为y(0) = 1。
在matlab中，可以通过以下代码来求解该微分方程：```定义微分方程的函数function dydt = myode(t, y)dydt = -y;调用ode45函数求解[t, y] = ode45(myode, [0, 5], 1);plot(t, y);```运行以上代码，即可得到微分方程的数值解，并通过绘图来展示解的变化。
六、总结迭龙格-库塔法是一种常用的数值解常微分方程的方法，它在matlab中有较为方便的调用方式。
通过ode45函数，可以快速求解常微分方程的数值解，并通过绘图来展示结果。
希望本篇文章对读者有所帮助，谢谢阅读。
七、应用场景和优势在实际应用中，迭龙格-库塔法广泛应用于各种科学和工程领域，如物理学、化学、生物学、经济学等。
四阶龙格库塔法求解微分方程组
四阶龙格库塔法求解微分方程组微分方程是数学中的一个重要分支，它描述了自然界中许多现象的发展规律。
在实际应用中，我们经常需要求解微分方程组，以得到系统的演化规律和性质。
本文将介绍一种常用的求解微分方程组的数值方法——四阶龙格库塔法。
一、微分方程组的数值解法微分方程组是形如下式的方程集合：begin{cases} frac{dy_1}{dx}=f_1(x,y_1,y_2,cdots,y_n) frac{dy_2}{dx}=f_2(x,y_1,y_2,cdots,y_n) cdotsfrac{dy_n}{dx}=f_n(x,y_1,y_2,cdots,y_n) end{cases} 其中，$y_1,y_2,cdots,y_n$是关于自变量$x$的未知函数，$f_1,f_2,cdots,f_n$是已知的函数。
求解微分方程组的数值方法主要有以下两种：1.欧拉法欧拉法是最简单的数值方法之一，它的基本思想是将微分方程组中的导数用差分代替，从而得到一个递推公式。
具体而言，欧拉法的递推公式为：$$y_{i+1}=y_i+hf(x_i,y_i)$$其中，$h$是步长，$x_i$和$y_i$分别是第$i$个点的自变量和因变量的值。
欧拉法的精度较低，容易出现数值误差，但是它的计算速度快，适用于一些简单的微分方程组。
2.龙格库塔法龙格库塔法是一种常用的高精度数值方法，它的基本思想是将微分方程组中的导数用一系列加权的差分代替，从而得到一个递推公式。
其中，四阶龙格库塔法是最为常用的一种。
具体而言，四阶龙格库塔法的递推公式为：begin{aligned} k_1&=hf(x_i,y_i)k_2&=hf(x_i+frac{h}{2},y_i+frac{k_1}{2})k_3&=hf(x_i+frac{h}{2},y_i+frac{k_2}{2})k_4&=hf(x_i+h,y_i+k_3)y_{i+1}&=y_i+frac{k_1+2k_2+2k_3+k_4}{6} end{aligned} 其中，$k_1,k_2,k_3,k_4$是加权的差分，$h$是步长，$x_i$和$y_i$分别是第$i$个点的自变量和因变量的值。
龙格-库塔法,求解常微分方程
隆格库塔法求解常微分方程摘要科学技术中常常需要求解常微分方程的定解问题，这里问题最简单的形式，是本章将着重考察的一阶方程的初值问题.虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些特殊类型的方程，实际问题中归结出来的微分方程主要靠数值解法求解.本文着重讨论了隆格库塔法求解一阶常微分方程的初值问题，采用了精度较高的经典的四阶隆格库塔法,然后通过对实例运用Matlab编程进行计算求解，为了体现计算结果的精确性和方法的优越性，再采用了欧拉法和预估较正法对实例进行计算求解作为比较.通过比较三种方法的计算精度，发现四阶经典龙格-库塔方法的误差最小，预估较正法其次，欧拉方法误差则比较大.最后通过选取不同的步长，研究了不同的步长对隆格库塔法求解常微分方程初值问题的计算精度的影响.总之，本文全面分析了隆格库塔法在求解常微分方程的应用，相比与其他的数值解法，隆格库塔法计算精度较高，收敛性较好，其中四阶的隆格库塔法的效率最高，精度也最高.关键词：四阶隆格库塔法；欧拉法；预估较正法；一阶常微分方程；MatlabRunge Kutta Method For Solving Ordinary Differential EquationsABSTRACTProblem solving ordinary differential equations are often needed in science andtechnology. the problem in the simplest form is the initial value problem of first order equations in this chapter ,which will be discussed. Although there are various analytical methods for solving ordinary differential equations, the analytical method can only be used to solve some special types of equations.differential equations can be summed up the actual problems whichThis paper discusses the initial value problem of Runge Kutta Barclays by solving a differential equation, using the four order Runge Kutta method with high accuracy.for instance through classic Matlab programming calculation, the superiority in order to accurately and reflect the calculation result, then the Euler method and the prediction correction method for instance by calculation through the calculation precision. The comparison of three kinds of methods, found that the error of four order Runge Kutta method of minimum, prediction correction method secondly, Euler method error is relatively large. Finally, by selecting different step, study the affect the calculation accuracy of different step of Runge Kutta method to solve initial value problems of ordinary differential equations.In short, this paper comprehensively analyzes the application of Runge Kutta method for solving ordinary differential equations, compared with the numerical solution of other, higher accuracy Runge Kutta method, good convergence, the Runge Kutta method of order four of the highest efficiency and its precision is the highest.Key words: Four order Runge Kutta method; Euler method; prediction correction method; first order ordinary differential equations; Matlab目录1 问题的提出 (1)1.1 问题背景............................................... . (1)1.2 问题的具体内容 (1)2 问题假设 (2)3 符号系统 (2)4 问题的分析 (3)4.1 欧拉格式 (3)4.2 预估较正法 (3)4.3 四阶隆格库塔法的格式 (4)5 模型的建立与求解 (4)5.1 隆格库塔法的基本原理 (4)5.1.1 Taylor级数 (4)5.1.2 隆格库塔法的基本思想 (4)5.1.3 四阶的隆格库塔法 (5)5.2 其他求解常微分方程边值问题算法的简介 (6)5.3 模型求解 (8)5.3.1 运用MATLAB软件对模型求解结果及析 (8)6 模型的评价 (16)7 课程设计的总结与体会 (16)参考文献 (17)附录 (18)一、问题的提出1.1 问题背景：科学技术中常常需要求常微分方程的定解问题，微分方程里最简单的方程形式莫过于一阶常微分方程的初值问题，即：0(,)()dy f x y a x b dx y a y ⎧=≤≤⎪⎨⎪=⎩ （1）其中a ，b 为常数.虽然求解此类微分方程有各种各样的解析方法，但解析方法只能用于求解一些特殊类型方程，实际问题中归结出来的微分方程主要靠数值解法求解.因为一阶常微分方程简单但又是求解其他方程的基础，所以发展了许多典型的解法.本文着重讨论一类高精度的单步法——隆格库塔法，并且运用四阶的隆格库塔格式来求解初值问题，并且通过实例运用四阶的隆格库塔格式来求解初值问题，同时与显式与隐式的Euler 格式求解出的结果进行精度比较.1.2 问题的具体内容实例一：在区间[0,1]上采用经典的四阶隆格库塔方法求解微分方程1(0)1dy y x dx y ⎧=-++⎪⎨⎪=⎩，其精确解为x y x e -=+，步长为0.5,然后用欧拉法，预估校正法分别求解，且将计算结果与精确解进行比较，对三个算法的收敛性的进行分析比较.实例二：在区间[0,1]上用经典的四阶龙格库塔方法求解初值问题 (0)1x dy xe y dx y -⎧=-⎪⎨⎪=⎩， 其精确解为21(2)2xx e -+，然后用欧拉法，预估校正法分别求解，且将计算结果与精确解进行比较，对三个算法的收敛性的进行分析比较.最后在区间[0,1]上分别取步长h=0.1;0.05时进行计算，并且探究选取不同的步长对计算结果精度的影响.二、问题假设2.1 假设数值方法本身的计算是准确的.2.2 假设选取的步长趋于0时计算的结果会收敛到微分方程的准确解.2.3 假设步长的增加不会导致舍入误差的严重积累.三、符号系统3.1 符号说明符号含义h选取的步长*K平均斜率p精度的阶数∆前后两次计算结果的偏差y第n个节点的实验值n()y x第n个节点的精确值nδ实验值与精确值的绝对误差四、问题的分析问题要求运用隆格库塔算法来求解一阶微分方程的初值问题，针对前面提出的实例，本文先用经典的四阶隆格库塔法来求解上面的微分方程，为了体现隆格库塔法的优越性，同时用欧拉法，预估校正法分别求解，且将计算结果与精确解进行比较，对三个算法的收敛性的进行分析比较.最后在区间[0,1]上分别取步长h=0.1;0.05时进行计算，分析在选取不同的步长时，求解结果的精度变化如何.下面是欧拉法，预估校正法以及经典的四阶隆格库塔法的计算公式.4.1欧拉格式（1）显式欧拉格式1(,)n n n n y y hf x y +=+ （2） 局部截断误差22211()()()()22n n n h h y x y y y x o h ξ++''''-=≈= （3） （2）隐式欧拉格式111(,)n n n n y y hf x y +++=+ （4）局部截断误差2211()()()2n n n h y x y y x o h ++''-≈-= (5) 4.2 预估校正法预估 1(,)n n n n y y hf x y +=+ (6)校正 111[(,)(,)]2n n n n n n h y y f x y f x y +++=++ (7) 统一格式 1[(,)(,(,))]2n n n n n n n n h y y f x y f x h y hf x y +=++++ (8) 平均化格式 11(,),(,),1().2p n n n c n n p n p c y y hf x y y y hf x y y y y ++⎧⎪=+⎪=+⎨⎪⎪=+⎩ (9)4.3 四阶龙格库塔方法的格式（经典格式）112341213243(22),6(,),(,),22(,),22(,).n n n n n n n n n n h y y K K K K K f x y h h K f x y K h h K f x y K K f x h y hK +⎧=++++⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩(10)五、模型的建立与模型求解5.1 隆格库塔法的基本原理隆格库塔法是一种高精度的单步法，这类方法与下述Taylor 级数法有着紧密的联系.5.1.1 Taylor 级数设初值问题 '00(,)()y f x y y x y ⎧=⎨=⎩有解，按泰勒展开，有2'''1()()()()....2n n n n h y x y x hy x y x +=+++； （11） 其中()y x 的各阶导数依据所给方程可以用函数f 来表达，下面引进函数序列(,)j f x y 来描叙求导过程，即(0)(1)'(0)''(1)(1)(1)'''(2),f f y f f y f f x x f f y f f x y ∂∂=≡=+≡∂∂∂∂=+≡∂∂ （12）(2)(2)()(1)j j j j f f y f f x y ---∂∂=+≡∂∂ （13） 根据上式，结果导出下面 Taylor 格式2'''()1...2!!pp n n nn n h h y y hy y y p +=++++ （14）其中一阶Taylor 格式为： '1n n n y y hy +=+提高Taylor 格式的阶数p 即可提高计算结果的精度，显然，p 阶Taylor 格式的局部截断误差为：11(1)1(1)!p p n n h y y y p ζ+++-+=+ （15） 因此它具有p 阶精度.5.1.2 隆格库塔方法的基本思想隆格库塔法实质就是间接地使用Taylor 级数法的一种方法,考察差商1()()n n y x y x h+- 根据微分中值定理，这有'1()()()n n n y x y x y x h h θ+-=+ （16） 利用所给方程 '(,)y f x y =1()()(,())n n n n y x y x hf x h y x h θθ+=+++ （17） 设 平均斜率*(,())n n K f x h y x h θθ=++，由此可见，只要对平均斜率一种算法，便相应地可以导出一种计算格式.再考察改进的Euler 格式，它可以改写成平均化的形式：1121211()2(,)(,)n n n n n n h y y K K K f x y K f x y hK ++⎧=++⎪⎪=⎨⎪=+⎪⎩（18） 这个处理过程启示我们，如果设法在1(,)n n x x +内多预测几个点的斜率值，然后将它们加权平均作为平均斜率，则有可能构造具有更高精度的计算格式，这就是隆格库塔法的基本思想.5.1.3 四阶的隆格库塔法为了方便起见，本文主要运用经典的隆格库塔算法-四阶隆格库塔格式.其格式如下：112341213243(22),6(,),(,),22(,),22(,).n n n n n n n n n n h y y K K K K K f x y h h K f x y K h h K f x y K K f x h y hK +⎧=++++⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩（19）下面为其具体的算法流程图：5.2 其他求解常微分边值问题算法的简介5.2.1欧拉数值算法（显式）微分方程里最简单的方程形式莫过于一阶常微分方程的初值问题，即：(,)()dyf x y a x b dx y a y ⎧=≤≤⎪⎨⎪=⎩(20) 其中a ，b 为常数.开始输入x 0,y 0,h,N x 1=x 0+hk 1=f(x 0,y 0),k 2=f(x 0+h/2,y 0+hk 1/2)k 3=f(x 0+h/2,y 0+hk 2/2),k 4=f(x 1,y 0+hk 3)y 1=y 0+h(k 1+2k 2+2k 3+k 4)/6n=1输出x 1,y 1n=N? 结束n=n+1x 0=x 1,y 0=y 1否是图5.1 龙格-库塔法流程图因为其简单但又是求解其他方程的基础，所以发展了许多典型的解法.所有算法中的f 就是代表上式中(,)f x y ，而y f 表示(,)f x y y ∂∂，x f 表示(,)f x y x∂∂. 简单欧拉法是一种单步递推算法.简单欧拉法的公式如下所示：1(,)n n n n y y hf x y +=+ (21)简单欧拉法的算法过程介绍如下:给出自变量x 的定义域[,]a b ，初始值0y 及步长h .对0,1,()/k b a h =-，计算1(,)k k k k y y hf x y +=+5.2.2欧拉数值算法（隐式）隐式欧拉法也叫退欧拉法，隐式欧拉法的公式如下所示：111(,)n n n n y y hf x y +++=+ (22)隐式欧拉法是一阶精度的方法，比它精度高的公式是：111[(,)(,)]2n n n n n n hy y f x y f x y +++=++ (23)隐式欧拉的算法过程介绍如下.给出自变量x 的定义域[,]a b ，初始值0y 及步长h .对0,1,()/k b a h =-，用牛顿法或其他方法求解方111(,)k k k k y y hf x y +++=+得出1k y +.5.2.3 欧拉预估-校正法改进欧拉法是一种二阶显式求解法，其计算公式如下所示：1[,(,)]22n n n n n n h h y y hf x y f x y +=+++11(,)[(,)(,)]2n n n n n n n n t y hf x y h y y f x y f x t ++=+⎧⎪⎨=++⎪⎩ (24)四阶龙格-库塔法有多种形式，除了改进的欧拉法外还有中点法.中点法计算公式为：1[,(,)]22n n n n n n h h y y hf x y f x y +=+++ (25)5.3 模型求解5.3.1运用MATLAB 软件对模型求解结果及分析用欧拉法、改进的欧拉格式、经典的四阶龙格库塔法来求解常微分方程的边值问题，并且比较其精度（具体的MATLAB 源程序见附录） 以下进行实例分析：实例一. 1(0)1dyy x dx y ⎧=-++⎪⎨⎪=⎩由题可知精确解为：x y x e -=+，当x=0时，y(x)=0.在这里取步长h 为0.1， 通过MATLAB 程序的计算，相应的结果如下：表5-1 步长为0.1时各方法的预测值与精确值的比较(精确到6位小数)初值 Euler 法 相对误差 预估校正法 相对误差 经典四阶库 相对误差精确值 0 -- -- -- -- -- -- 1.00000 0.1 1.00910 0.00424 1.00500 0.00016 1.00484 0.00000 1.00484 0.2 1.02646 0.00759 1.01903 0.00029 1.01873 0.00000 1.01873 0.3 1.05134 0.01011 1.04122 0.00038 1.04082 0.00000 1.04082 0.4 1.08304 0.01189 1.07080 0.00045 1.07032 0.00000 1.07032 0.5 1.12095 0.01303 1.10708 0.00049 1.10653 0.00000 1.10653 0.6 1.16451 0.01366 1.14940 0.00052 1.14881 0.00000 1.14881 0.7 1.21319 0.01388 1.19721 0.00052 1.19659 0.00000 1.19659 0.8 1.26655 0.01378 1.24998 0.00052 1.24933 0.00000 1.24933 0.9 1.32414 0.01344 1.30723 0.00050 1.30657 0.00000 1.30657 1.01.38558 0.01294 1.36854 0.00048 1.36788 0.000001.36788步长为0.1时的精确值与预测值的比较精确值欧拉法改进欧拉格式四阶龙格库塔轴Y00.10.20.30.40.50.60.70.80.91 1.1 1.2 1.3 1.4X轴图5.2 步长为0.1时各方法的预测值与精确值的比较原函数图像轴YX轴图5.3 步长为0.1时原函数图像在这里取步长h为0.05，通过MATLAB程序的计算，相应的结果如下：表5-2 h=0.05时三个方法与精确值的真值表步长Euler法相对误差预估校正法相对误差经典四阶库相对误差精确值0 -- -- -- -- -- --1.00000 0.05 1.00250 0.00911 1.00125 0.01035 1.00123 0.01037 1.00484 0.10 1.00738 0.01711 1.00488 0.01954 1.00484 0.01958 1.01873 0.15 1.01451 0.02405 1.01076 0.02765 1.01071 0.02770 1.04082 0.20 1.02378 0.03001 1.01880 0.03473 1.01873 0.03479 1.07032 0.25 1.03509 0.03507 1.02889 0.04085 1.02880 0.04093 1.10653 0.30 1.04834 0.03930 1.04092 0.04610 1.04082 0.04619 1.14881 0.35 1.06342 0.04277 1.05480 0.05053 1.05469 0.05063 1.19659 0.40 1.08025 0.04555 1.07044 0.05422 1.07032 0.05432 1.24933 0.45 1.09874 0.04772 1.08775 0.05724 1.08763 0.05734 1.30657 0.50 1.11880 0.04933 1.10666 0.05964 1.10653 0.05975 1.36788 0.55 1.14036 0.05045 1.12709 0.06150 1.12695 0.06161 1.00484 0.60 1.16334 0.05113 1.14895 0.06286 1.14881 0.06298 1.01873 0.65 1.18768 0.05143 1.17219 0.06379 1.17205 0.06391 1.04082 0.70 1.21329 0.05139 1.19674 0.06433 1.19659 0.06445 1.07032 0.75 1.24013 0.05106 1.22252 0.06453 1.22237 0.06465 1.10653 0.80 1.26812 0.05048 1.24949 0.06443 1.24933 0.06455 1.14881 0.85 1.29722 0.04969 1.27757 0.06408 1.27742 0.06419 1.19659 0.90 1.32735 0.04871 1.30673 0.06349 1.30657 0.06361 1.249330.95 1.35849 0.04759 1.33690 0.06272 1.33674 0.06283 1.306571.00 1.39056 0.04634 1.36804 0.06178 1.36788 0.06189 1.36788欧拉格式 改进欧拉格式四阶龙格库塔精确值图5.4 步长为0.05时各方法的预测值与精确值的比较11.051.11.151.21.251.31.351.4X 轴Y 轴原函数图像图5.5 步长为0.05时原函数图像实例2 (0)1xdy xe ydx y -⎧=-⎪⎨⎪=⎩由题可知精确解为：21(2)2x x e -+当x=0时，y(x)=0.在这里取步长h为0.1，通过MATLAB程序的计算，相应的结果如下：步长Euler法相对误差预估校正法相对误差经典四阶库相对误差精确值0.1 0.9000 0.0318 0.9096 0.0214 0.9094 0.0216 0.9295 0.2 0.8192 0.0611 0.8359 0.0420 0.8356 0.0424 0.8726 0.3 0.7544 0.0876 0.7761 0.0614 0.7757 0.0619 0.8268 0.4 0.7027 0.1109 0.7277 0.0793 0.7272 0.0799 0.7903 0.5 0.6617 0.1310 0.6886 0.0956 0.6881 0.0963 0.7615 0.6 0.6294 0.1480 0.6572 0.1103 0.6567 0.1110 0.7387 0.7 0.6040 0.1621 0.6320 0.1234 0.6315 0.1241 0.7209 0.8 0.5841 0.1737 0.6116 0.1349 0.6111 0.1355 0.70690.9 0.5686 0.1830 0.5951 0.1450 0.5946 0.1456 0.69591.0 0.5563 0.1905 0.5815 0.1538 0.5811 0.1544 0.6872原函数图像轴Y00.10.20.30.40.50.60.70.80.91X轴图5.6 步长为0.1时原函数图像各方法的预测值与精确值的比较欧拉格式改进的格式四阶龙格库塔精确值轴Y0.10.20.30.40.50.60.70.80.91X轴图5.7 步长为0.1时各方法的预测值与精确值的比较在这里取步长h为0.05，通过MATLAB程序的计算，相应的结果如下：表5-4 步长为0.1时各方法的预测值与精确值的比较(精确到5位小数)步长Euler法相对误差预估校正法相对误差经典四阶库相对误差精确值0.050.95000 0.01342 0.95245 0.01088 0.95242 0.01090 0.96292 0.100.90490 0.02650 0.90947 0.02158 0.90942 0.02163 0.92953 0.150.86428 0.03916 0.87067 0.03206 0.87061 0.03213 0.89951 0.200.82774 0.05137 0.83567 0.04228 0.83559 0.04237 0.87256 0.250.79491 0.06307 0.80414 0.05219 0.80405 0.05230 0.84842 0.300.76545 0.07423 0.77576 0.06176 0.77566 0.06188 0.82682 0.350.73904 0.08482 0.75023 0.07096 0.75013 0.07110 0.80754 0.40 0.71541 0.09482 0.72730 0.07977 0.72719 0.07991 0.79035 0.45 0.69427 0.10422 0.70672 0.08817 0.70660 0.08832 0.77505 0.50 0.67539 0.11302 0.68825 0.09615 0.68813 0.09630 0.76146 0.55 0.65855 0.12123 0.67168 0.10370 0.67156 0.10386 0.74940 0.60 0.64352 0.12886 0.65683 0.11084 0.65671 0.11100 0.73871 0.65 0.63012 0.13593 0.64351 0.11756 0.64340 0.11773 0.72925 0.70 0.61819 0.14245 0.63157 0.12388 0.63145 0.12405 0.72087 0.75 0.60754 0.14847 0.62085 0.12981 0.62073 0.12998 0.71347 0.80 0.59805 0.15400 0.61121 0.13537 0.61110 0.13553 0.70691 0.85 0.58957 0.15908 0.60254 0.14058 0.60243 0.14073 0.70109 0.90 0.58198 0.16374 0.59470 0.14545 0.59460 0.14560 0.695930.95 0.57517 0.16802 0.58761 0.15001 0.58751 0.15016 0.691321.00 0.56904 0.17194 0.58117 0.15429 0.58107 0.15443 0.687190.10.20.30.40.50.60.70.80.910.550.60.650.70.750.80.850.90.951X 轴Y 轴各方法下的预测值与精确值的比较欧拉格式改进欧拉格式四阶龙格库塔精确值图5.8 步长为0.05时各方法的预测值与精确值的比较六、模型的评价本文着重讨论了4阶的隆格库塔法来求解微分方程，并且通过两个实例验证了隆格库塔法在求解初值问题的优越性.从上面的实例比较可知，在计算精度上，四阶经典龙格-库塔方法的误差最小，改进欧拉方法其次，欧拉方法误差则比较大，所以四阶经典龙格-库塔方法得到最佳的精度.而在计算量上面，相应地，很明显的四阶经典龙格-库塔方法也是最大，改进欧拉方法其次，欧拉方法计算量最小.这样的结果，说明了运用以上三种方法时，其计算量的多少与精度的大小成正比.我们在实际运用与操作中，可以根据实际情况，选择这3种方法中的其中一种最适合的，追求精度的话，可以使用四阶经典龙格-库塔方法；而改进的欧拉方法，在精度上和计算量上都表现得很出色，能够满足一般情况；而欧拉方法更主要的是适用于对y的估计上，而精度则有所欠缺，以上各方法的选择，都取决于具体的情况.七、课程设计的总结与体会本文着重采用隆格库塔法运用MATLAB编程来求解微分方程，相比于欧拉法以及预估校正法，隆格库塔法在提高近似值解的精度上是非常起作用的，而且又具有计算量不大、算法组织容易.其次，每一次的课程设计总是让我学到了更多的知识，不论是C++、SPASS还是MATLAB软件，这些都让我学到了如何解决实际问题的好工具，通过这些工具，是自己能够得到突破和成长.以下是我完成此次课程设计的几点体会：（1）必须学好基础知识，在做的过程中，发现自己有很多东西都不懂，要博学必须从一点一点做起.以往训练得少只是把握的不牢靠.所以做起来感到有点吃力.所以，无论什么学科，一定要打好基础.（2）程序设计要靠多练，多见识，那样形成一种编程思维，我想对我是有很大好处的.尤其像我这种平时学得不扎实的人.（3）做事情要有恒心，遇到困难不要怕，坚决去做.如果做出来了，固然高兴，如果没有做出来也没关系，自己努力了对得起自己就好.同时，把它看做是对自己的锻炼. （4）做程序特别是做大程序是很有趣的.虽然有的问题很难，要花很多时间很多精力，但是那种解决了一个问题时的喜悦足以把付出的辛苦补偿回来.得到一种心里的慰藉.参考文献[1] 李庆杨，王能超，易大义编．数值分析（第四版）[M]：华中科技大学出版社，2006.[2] 姜启源，谢金星，叶俊编．数学模型（第三版）[M]．北京：高等教育出版社，2005[3] 刘琼荪，数学实验[M]，高等教育出版社，2004[4] 王建伟，MATLAB7.X程序设计[M]，中国水利水电出版社，2007[5] 王高雄，周之铭等编.常微分方程（第三版） [M]:高等教育出版社，2006[6]何坚勇编著. 运筹学基础（第二版）[M]. 北京：清华大学出版社，2008.附录附录一：显示欧拉法matlab程序%欧拉法clear allclcx=[];y=[];y1=[];h=0.1;x=0:h:1;n=length(x);for i=1:ny(i)=f1(x(i));endfigure(1)plot(x,y,'g-');hold ony1(1)=1;for j=2:ny1(j)=y1(j-1)+h*f(x(j-1),y1(j-1)); endY=[x;y1];fid=fopen('data.txt','wt');fprintf(fid,'%6.2f %12.4f\n',Y);fclose(fid);plot(x,y1,'r-');figure(2)DT=abs(y-y1);plot(x,DT)%1.建立导数函数文件function z=f(x,y)z=y-2*x/y;%2.建立原函数文件function z1=f1(x)z1=(2*x+1)^(1/2);迭代n次的后退的欧拉格式matlab程序%1.建立导数函数文件function z=f(x,y)z=y-2*x/y;%2.建立原函数文件function z1=f1(x)z1=(2*x+1)^(1/2);%迭代n次的后退的欧拉格式clear allclcN=input('请输入迭代次数:');x=[];y=[];y1=[];h=0.1;x=0:0.1:1;n=length(x);for i=1:ny(i)=f1(x(i));endfigure(1)plot(x,y,'g-');for j=2:ny1(j)=y1(j-1)+h*f(x(j-1),y1(j-1)); T=y1(j);for k=1:NT=y1(j-1)+h*f(x(j),T);endy1(j)=T;endY=[x;y1];fid=fopen('data.txt','wt');fprintf(fid,'%6.2f %12.4f\n',Y); fclose(fid);plot(x,y1,'r-');figure(2)DT=abs(y-y1);plot(x,DT)附录二：预估校正法matlab程序%1.建立导数函数文件function z=f(x,y)z=y-2*x/y;%2.建立原函数文件function z1=f1(x)z1=(2*x+1)^(1/2);%预估校正法%欧拉法clear allclcx=[];y1(1)=1;y2=[];y2(1)=1;h=0.1;x=0:0.1:1;n=length(x);for i=1:ny(i)=f1(x(i));endfigure(1)plot(x,y,'g-');hold onfor j=2:ny1(j)=y2(j-1)+h*f(x(j-1),y2(j-1));y2(j)=y2(j-1)+(h/2).*[f(x(j-1),y2(j-1))+f(x(j),y1(j))]; endY=[x;y2];fid=fopen('data.txt','wt');fprintf(fid,'%6.2f %12.4f\n',Y);fclose(fid);plot(x,y2,'r-');figure(2)DT=y-y2;plot(x,DT)附录三：四阶龙格库塔法matlab程序%四阶龙格库塔法clear allclcn=length(x);y=[];y1=[];y1(1)=1;for i=1:ny(i)=f1(x(i));endfigure(1)plot(x,y,'g-');hold onfor j=2:nK1=f(x(j-1),y1(j-1));K2=f(x(j-1)+h/2,y1(j-1)+h/2*K1);K3=f(x(j-1)+h/2,y1(j-1)+h/2*K2);K4=f(x(j-1)+h,y1(j-1)+h*K3);y1(j)=y1(j-1)+(h/6)*(K1+2*K2+2*K3+K4); endY=[x;y1];fid=fopen('data1.txt','wt');fprintf(fid,'%6.2f %12.4f\n',Y);fclose(fid);plot(x,y1,'r-');figure(2)DT=abs(y-y1);plot(x,DT)。
龙格库塔
龙格-库塔(Runge-kutta)算法在MATLAB中的用法微分方程ode23:使用二阶龙格-库塔法求微分方程，调用格式为：[t，y]=ode23(‘fname’,tspan,y0)ode45:使用四阶龙格-库塔法求微分方程，调用格式为：[t，y]=ode45(‘fname’,tspan,y0)其中fname为由M函数定义的线性或者非线性微分方程的句柄函数名。
tspan形式为[t0,tf]，表示求解区间。
y0是初始状态列向量。
t和y分别给出时间向量和相应的状态向量。
二阶龙格-库塔法(ode23)：下面式2为Euler(欧拉法)增量，即一步起始端斜率，式3为一步终点端斜率。
所以式1仿真计算的增量实际上是取两点斜率的平均斜率来计算的，其精度高于Euler算法。
四阶龙格-库塔法(ode45)：计算原理为预报-校正法，预报值采用Euler算出，下式又作了3次校正，因此计算精度远高于Euler算法。
e.g.二阶非线性系统的微分方程：x″ + 0.5*x′+ 2*x + x^2 = 0求系统在初始条件为x(0)=1，x′(0)=0的数值解。
建立M函数：function xdot=odetest(t,x)%龙格-库塔算法测试%二阶非线性系统微分方程xdotdot + 0.5*xdot + 2*x + x^2 = 0 %求系统在初始条件为x(0)=1,xdot(0)=0的数值解xdot=zeros(2,1);xdot(1)=x(2);xdot(2)=-0.5*x(2)-2*x(1)-x(1)^2;tspam=[0 20];x0=[0 1];[t,x]=ode23(@odetest,tspan,x0);plot(t,x)微分方程数值解曲线如下图：。
runge-kutta方法
runge-kutta方法Runge-Kutta方法是求解常微分方程组的一种常用数值方法，它是由德国数学家卡尔·龙格和马丁·库塔（Martin Kutta）发明的。
该方法是一种高阶精度的方法，可以求解一阶或高阶的常微分方程组。
$$ y_{n+1} = y_n + \sum_{i=1}^s b_i k_i $$$y_n$表示第$n$个步骤中的输出，$y_{n+1}$是下一个步骤的输出。
$k_i$表示第$i$个中间变量，$b_i$是权值，它们是根据所采用的方法计算得出的。
我们将在下文中详细介绍这些变量和权值的计算方法。
在数值计算中，我们通常采用自适应步长的方法，这需要我们给出一个误差容限值$\epsilon$，根据误差限定步长大小。
如果当前步长不满足误差要求，我们会将步长缩小再进行计算，或者选择更高阶的方法进行计算，以提高精度。
接下来我们将详细介绍如何使用Runge-Kutta方法求解常微分方程组。
1. Runge-Kutta方法的中间变量$k_i$的计算我们需要计算中间变量$k_i$的值。
它们的计算公式如下：$$ k_1 = hf(t_n, y_n) $$$$ k_2 = hf(t_n + c_2 h, y_n + a_{21}k_1) $$$$ ... $$$$ k_s = hf(t_n + c_s h, y_n + \sum_{j=1}^{s-1}a_{sj}k_j) $$$h$表示步长大小，$t_n$和$y_n$分别表示当前时间和状态，$c_i$和$a_{ij}$是与所采用的方法相关的常数。
$$ c_2 = \frac{1}{2} , a_{21} = \frac{1}{2}$$- 改进的Runge-Kutta方法（RK4）：根据具体选择的方法，我们可以计算出$k_i$的值。
$b_{ij}$也是与所采用的方法相关的常数，对于不同的方法而言，它们的计算公式不同。
根据计算出的中间变量$k_i$和权值$b_i$，我们可以进行Runge-Kutta方法的迭代，求解常微分方程组。
matlab用四阶龙格库塔函数求解微分方程组
一、介绍Matlab作为一种强大的科学计算软件，提供了众多函数和工具来解决微分方程组。
其中，四阶龙格库塔函数是一种常用的数值方法，用于求解常微分方程组。
本文将介绍如何使用Matlab中的四阶龙格库塔函数来求解微分方程组，并对该方法的原理和实现进行详细说明。
二、四阶龙格库塔方法四阶龙格库塔方法是一种常用的数值方法，用于求解常微分方程组。
它是一种显式的Runge-Kutta方法，通过逐步逼近微分方程的解，在每一步使用多个中间值来计算下一步的解。
该方法通过四个中间值来计算下一步的状态，并且具有较高的精度和稳定性。
三、在Matlab中使用四阶龙格库塔方法求解微分方程组在Matlab中，可以使用ode45函数来调用四阶龙格库塔方法来解决微分方程组的问题。
ode45函数是Matlab提供的用于求解常微分方程组的函数，可以通过指定微分方程组以及初值条件来调用四阶龙格库塔方法来进行求解。
1. 定义微分方程组我们需要定义要求解的微分方程组。
可以使用Matlab中的匿名函数来定义微分方程组，例如：```matlabf = (t, y) [y(2); -sin(y(1))];```其中，f是一个匿名函数，用于表示微分方程组。
在这个例子中，微分方程组是y' = y2, y2' = -sin(y1)。
2. 指定初值条件和求解区间接下来，我们需要指定微分方程组的初值条件和求解区间。
初值条件可以通过指定一个初始时刻的状态向量来完成，例如：```matlabtspan = [0, 10];y0 = [0, 1];```其中，tspan表示求解区间，y0表示初值条件。
3. 调用ode45函数进行求解我们可以通过调用ode45函数来求解微分方程组的数值解。
具体的调用方式如下：```matlab[t, y] = ode45(f, tspan, y0);```其中，t和y分别表示求解的时间点和对应的状态值。
四、示例下面我们通过一个具体的例子来演示如何使用Matlab中的四阶龙格库塔方法来求解微分方程组。
matlab用四阶龙格库塔法解二阶常微分方程
matlab用四阶龙格库塔法解二阶常微分方程在数学和工程中，常微分方程是描述自然界中各种物理现象和过程的常见数学模型。
常微分方程通常包含未知函数及其导数的方程。
在求解常微分方程时，人们通常使用数值方法来近似求解，其中一种常见的方法是使用龙格-库塔方法。
四阶龙格-库塔方法是一种常用的数值方法，用于求解二阶常微分方程。
它是通过在每个步骤中估计未知函数的导数来数值求解方程。
它的精度比较高，通常被认为是最常用的龙格-库塔方法。
对于一个二阶常微分方程：y''(t) = f(t, y(t), y'(t))其中y(t)是未知函数，f(t, y(t), y'(t))是已知的函数。
我们可以将这个二阶微分方程转化为两个一阶微分方程：y'(t) = u(t)u'(t) = f(t, y(t), u(t))其中u(t)是y'(t)的一个替代函数。
为了使用四阶龙格-库塔方法求解这个方程，我们需要将时间区间[t0, tn]分割为等距的n个子区间，每个子区间的长度为h = (tn - t0)/n。
我们从初始点t0开始，通过多个步骤来逼近解直到tn。
每一步我们使用以下递推公式来估计y(t)和u(t)的值：k1 = h * u(t)l1 = h * f(t, y(t), u(t))k2 = h * (u(t) + l1/2)l2 = h * f(t + h/2, y(t) + k1/2, u(t) + l1/2)k3 = h * (u(t) + l2/2)l3 = h * f(t + h/2, y(t) + k2/2, u(t) + l2/2)k4 = h * (u(t) + l3)l4 = h * f(t + h, y(t) + k3, u(t) + l3)然后我们可以使用以下公式来更新y(t)和u(t)的值：y(t + h) = y(t) + (k1 + 2k2 + 2k3 + k4)/6u(t + h) = u(t) + (l1 + 2l2 + 2l3 + l4)/6这个过程重复n次，直到达到结束时间tn。
求解随机微分方程的三级半隐式随机龙格库塔方法
求解随机微分方程的三级半隐式随机龙格库塔方法随机微分方程是具有随机项的微分方程，它在许多领域的研究中发挥着重要的作用。
随机微分方程的数值解法是研究中的一个重要问题，其中随机龙格库塔方法是常用的一种数值解法之一、本文将介绍随机微分方程的一种三级半隐式随机龙格库塔方法。
首先，我们考虑如下形式的随机微分方程：$$dX(t) = a(t,X(t))dt + b(t,X(t))dW(t)$$其中，$X(t)$是未知的随机过程，$a(t,X(t))$和$b(t,X(t))$是已知函数，$W(t)$是一个标准布朗运动。
我们的目标是求解方程在给定的时间间隔$[0,T]$内的数值解。
为了进行时间离散化，我们将时间间隔[0, T]分成N个小时间步长$\Delta t = \frac{T}{N}$。
令$t_i = i\Delta t$，$i = 0,1,2,...,N$，我们可以将方程改写为：$$X(t_{i+1}) = X(t_i) + a(t_i,X(t_i))\Delta t +b(t_i,X(t_i))\Delta W_i$$其中，$\Delta W_i = W(t_{i+1})-W(t_i)$是布朗运动在时间步长$\Delta t$内的增量。
注意到在上式中，$X(t_{i+1})$是未知的，我们需要进行反复迭代求解。
为了简化计算，我们引入半隐式随机龙格库塔方法。
半隐式随机龙格库塔方法将一阶随机微分方程以二阶精度数值求解，其中随机项以前一时间步长$t_i$的值来近似。
在本文中，我们将介绍一种三级半隐式随机龙格库塔方法，采用其中一种方式来估计方程的解。
首先，我们将时间$t$的导数项$a(t,X(t))$以及随机项$b(t,X(t))$在时间步$t_i$进行泰勒展开：$$a(t,X(t)) = a(t_i,X(t_i)) + \frac{\partiala(t,X(t))}{\partial t}，_{t_i} (t_{i+1} - t_i) + \frac{\partiala(t,X(t))}{\partial X}，_{t_i} (X(t_{i+1}) - X(t_i)) + O(\Deltat^2)$$$$b(t,X(t)) = b(t_i,X(t_i)) + \frac{\partialb(t,X(t))}{\partial t}，_{t_i} (t_{i+1} - t_i) + \frac{\partialb(t,X(t))}{\partial X}，_{t_i} (X(t_{i+1}) - X(t_i)) + O(\Deltat^2)$$将上述展开式代入原方程，我们可以得到：$$X(t_{i+1}) = X(t_{i}) + (a(t_i,X(t_i)) + \frac{\partiala(t,X(t))}{\partial X}，_{t_i} (X(t_{i+1}) - X(t_i)))\Delta t + (b(t_i,X(t_i)) + \frac{\partial b(t,X(t))}{\partial X}，_{t_i} (X(t_{i+1}) - X(t_i)))\Delta W_i$$接下来，我们采用不同方式来估计方程的解。
随机常微分方程的龙格库塔解法
随机常微分方程的龙格库塔解法
龙格库塔解法是一种用于解决随机常微分方程的常用方法。
它是一种把随机微分方程分解成非线性方程组的近似解法，通常可以用来解决系统的非线性演化方程，其中每个方程都有随机性。
它是一种被广泛应用于自然科学和工程领域的数值解法，它可以用来求解随机性较强的系统和不稳定性较强的系统的动力学行为。
龙格库塔解法的基本思想是将随机常微分方程拆分成一系列的近似子问题，从而使得系统的动力学行为可以精确的描述。
它通过将方程中的随机变量进行离散化，将复杂的随机微分方程转换为一系列的近似子问题，然后通过解决这些子问题来求解原始随机微分方程。
龙格库塔解法有一定的计算复杂度，但是它具有较高的精度，能够有效地描述系统的动力学行为。
它还具有较好的可扩展性，能够有效地解决复杂问题。
总之，龙格库塔解法是一种用于解决随机常微分方程的有效方法，它具有精度高，可扩展性好，计算复杂度低，并且广泛应用于许多研究领域的特点。
因此，它被认为是一种有效的数值解法，可以用来模拟复杂的系统，并得到准确的结果。
欧拉法、梯形法和龙格_库塔解微分方程
欧拉法、梯形法和龙格-库塔一、解方程：= 8x(2-y)y(0)=1二、算出方程的解析解为：y= 2 -三、实验原理：1.欧拉法原理：将区间[a,b]分成n段.那么方程在第x i点有y'(x i) = f(x i,y(x i)).再用向前差商近似代替导数则为：.在这里.h是步长.即相邻两个结点间的距离。
因此可以根据xi点和yi点的数值计算出y i+1来：.i=0,1,2,n这就是欧拉格式.若初值y i+ 1是已知的.则可依据上式逐步算出数值解y1,y2……..yn2.梯形法原理：将向前欧拉公式中的导数f(xi,yi)改为微元两端f(xi,yi)和f(xi+1,yi+1)的平均.即梯形公式。
3.龙格-库塔方法的基本思想：在区间[xn,xn+1]内多取几个点.将他们的斜率加权平均.作为导数的近似。
令初值问题表述如下。
则.对于该问题的RK4由如下方程给出：其中这样.下一个值(y n+1)由现在的值(y n)加上时间间隔(h)和一个估算的斜率的乘积决定。
该斜率是以下斜率的加权平均：k1是时间段开始时的斜率；k2是时间段中点的斜率.通过欧拉法采用斜率k1来决定y在点tn + h/2的值；k3也是中点的斜率.但是这次采用斜率k2决定y值；k4是时间段终点的斜率.其y值用k3决定。
当四个斜率取平均时.中点的斜率有更大的权值：RK4法是四阶方法.也就是说每步的误差是h5阶.而总积累误差为h4阶。
四、欧拉法、梯形法和龙格-库塔的实现代码：h=0.1;x=0:h:1;y1=zeros(size(x));y1(1)=1;y2=zeros(size(x));y2(1)=1;y3=zeros(size(x));y3(1)=1;for i1=2:length(x)y1(i1) = y1(i1-1) + h*8*x(i1-1)*(2-y1(i1-1));%欧拉法 m1= 8*x(i1-1)*(2-y2(i1-1));%梯形法m2 =8*x(i1)*(2-y2(i1-1)+h*m1);y2(i1)=y2(i1-1)+h*(m1+m2)/2;%梯形法公式k1=8*x(i1-1)*(2-y2(i1-1)); %龙格-库塔k2=8*(x(i1-1)+h/2)*(2-(y2(i1-1)+h*k1/2));k3=8*(x(i1-1)+h/2)*(2-(y2(i1-1)+h*k2/2));k4=8*(x(i1-1)+h)*(2-(y2(i1-1)+h*k3));y3(i1)=y3(i1-1)+(k1+2*k2+2*k3+k4)*h/6; %龙格-库塔公式endy4=2-exp(-4*(x.^2));%解析解plot(x,y1,x,y2,x,y3,x,y4)%解析解与数值解图像legend('y1','y2','y3','y4')plot(x,y4-y1,x,y4-y2,x,y4-y3)% 解析解与数值解误差图像legend('y4-y1','y4-y2','y4-y3')五、图像：1.解析解y4与各个数值解y1,y2,y3的图像:(1)当步长h=0.1时：（2）当步长h=0.05时：（3）当步长h=0.01时：结论：三个图中.方程的解析解都是y4。
龙格-库塔求解微分方程
MATLAB 在数值分析中的应用摘 要：数值分析是研究数学问题数值解及其理论的一个数学分支,涉及面很广.自计算机在其中应用以后,使其应用价值更为广泛,解决实际问题更为有效.以MA TLAB 为平台的数值分析方法与图形可视化编程在解决一些复杂问题时,大大降低了难度,减少了工作量.本文以经典龙格-库塔算法为核心,介绍用数值方法解常微分方程(组),分析自制系统的稳定性问题.如范德波尔方程.以及混沌现象中奇怪吸引子的模拟实现关键词： MA TLAB; 经典龙格-库塔算法; 解常微分方程(组); 范德波尔方程; 奇怪吸引子;1 MATLAB 平台1.1 矩阵MATLAB 是矩阵实验室（Matrix Laboratory ）的简称, 美国MathWorks 公司出品的商业数学软件.以下介绍在处理方程组时的特殊技法: A=[a b c; d e f; g h i]得到矩阵a b c d e f g h i也可直接输入方程组:function y=equition(t,x)y=[q(2)+q(1)/sqrt(q(1)^2+q(2)^2)*(1-q(1)^2-q(2)^2); % 以q 向量保存方程组;-q(1)+q(2)/sqrt(q(1)^2+q(2)^2)*(1-q(1)^2-q(2)^2)] % q(1)代表应变量x, q(2) 代表应变量y得到1.2 语言MATLAB 编程语言与 C 语言 类似, 格式上略有不同. 里不做详述, 程序中会做必要的标注.但在处理矩阵式十分方便, C++是以数组的形式存放变量,对某个元素都要准确指出其位置.22222222()]()]x y x y x y y x x y x y ⎧=+-+⎪+⎪⎨⎪=-+-+⎪+⎩C++是以数组的形式存放变量,对某个元素都要准确指出其位置,处理循环时算法复杂,即使使用指针也很容易出错.2经典龙格-库塔算法公式(15张)龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。
数值计算：四阶龙格-库塔法for二阶微分方程
数值计算：四阶龙格-库塔法for⼆阶微分⽅程引⾔考虑存在以下⼆阶偏微分⽅程f2⋅¨X(t)+f1⋅˙X(t)+f0⋅X(t)=F(t)如何使⽤四阶龙格-库塔法求解该微分⽅程？⼀阶微分⽅程的解法⾸先回顾下对于⼀阶微分⽅程的解法，现在有以下⼀阶常系数⾮齐次微分⽅程f1⋅˙X(t)+f0⋅X(t)=F(t)⽅程可以写作˙X(t)=F(t)−f0⋅X(t)f1X n+1−X ndt=F(t)−f0⋅X(t)f1X n+1=X n+dt⋅F(t)−f0⋅X nf1=X n+dt⋅f(t,X n)其中dt为时间步长。
四阶龙格-库塔法公式如下：X n+1=X n+dt6k1+2k2+2k3+k4k1=f t,X nk2=f t+dt2,Xn+dt2⋅k1k3=f t+dt2,Xn+dt2⋅k2k4=f t+dt,X n+dt⋅k3利⽤以上公式即可显式推进求解⼆阶微分⽅程的解法现在考虑⼆阶常系数⾮齐次微分⽅程f2⋅¨X(t)+f1⋅˙X(t)+f0⋅X(t)=F(t)⽅程中出现⼆次导数项¨X(t)，难以处理，所以考虑将⼆次导数项¨X(t)转化为⼀次导数项，⾄此，我们引⼊a(t)=X(t)b(t)=˙a(t)=˙X(t)原⽅程可以写成⼀阶偏微分⽅程组{()()()()(){f 2⋅˙b (t )+f 1⋅b (t )+f 0⋅a (t )=F (t )b (t )=˙a(t )仿照以上⼀阶微分⽅程的RK4解法˙a(t )=b (t )˙b (t )=F (t )−(f 1⋅b (t )+f 0⋅a (t ))f 2˙a n +1=a n +dt ⋅b (t )=a n +dt ⋅f (t ,a n ,b n )˙bn +1=b n +dt ⋅F (t )−(f 1⋅b (t )+f 0⋅a (t ))f 2=b n +dt ⋅g (t ,a n ,b n )˙a n +1=a n +dt ⋅f (t ,a n ,b n )˙bn +1=b n +dt ⋅g (t ,a n ,b n )化简⾄此，便可以仿照之前的RK4⽅法进⾏求解，具体公式为a n +1b n +1=a n +dt6⋅(k 1a +2k 2a +2k 3a +k 4a )b n +dt6⋅(k 1b +2k 2b +2k 3b +k 4b )其中，各个系数求解公式为：k 1a k 1b=f (t ,a n ,b n )g (t ,a n ,b n )k 2a k 2b =f (t +dt2,a n +dt2⋅k 1a ,b n +dt2⋅k 1b )g (t +dt2,a n +dt2⋅k 1a ,b n +dt2⋅k 1b )k 3a k 3b=f (t +dt2,a n +dt2⋅k 2a ,b n +dt2⋅k 2b )g (t +dt2,a n +dt2⋅k 2a ,b n +dt2⋅k 2b )k 4a k 4b=f (t +dt ,a n +dt ⋅k 3a ,b n +dt ⋅k 3b )g (t +dt ,a n +dt ⋅k 3a ,b n +dt ⋅k 3b )⾄此已经完成使⽤四阶龙格-库塔法⼆阶微分⽅程求解过程的梳理实例求解，以下偏微分⽅程：\begin{align} 1 \cdot \ddot{X(t)}+0 \cdot \dot{X(t)} +0 \cdot {X(t)} =cos(t)\\ \dot{X(0)}=0\\ {X(0)} =0 \end{align}理论解：\begin{align} {X(t)} =-cos(t)+1 \end{align}#include"RK_ODE.h"#include<iostream>using namespace std;#define M 1{{{{[][]{[][][][][][][][]MatrixXd F2(double t) {MatrixXd res = MatrixXd::Identity(M, M);return res;}MatrixXd F1(double t) {MatrixXd res = MatrixXd::Zero(M, M);return res;}MatrixXd F0(double t) {MatrixXd res = MatrixXd::Zero(M, M);return res;}MatrixXd F(double t) {MatrixXd res = MatrixXd::Ones(M, 1);res(0, 0) = cos(t);return res;}int main() {RK_ODE ode(M, 0.1);ode.X.fill(0.);ode.dX.fill(0.);for (int i = 0; i < 1000; i++) {ode.Solve(F2, F1, F0, F);ode.WriteMotionAndForce("results1.txt", 0); }//ode.SaveSnapshoot("snapshoot.json");cout << ode._Mode << endl;}Processing math: 88%。
四阶龙格库塔法解微分方程
四阶龙格库塔法解微分方程一、四阶龙格库塔法解一阶微分方程①一阶微分方程：cos y t ，初始值y(0)=0，求解区间[0 10]。
MATLAB 程序：%%%%%%%%%%% 四阶龙哥库塔法解一阶微分方程%%%%%%%%%%% y'=cost%%%%%%%%%%% y(0)=0, 0≤t ≤10，h=0.01%%%%%%%%%%% y=sinth=0.01;hf=10;t=0:h:hf;y=zeros(1,length(t));y(1)=0;F=@(t,y)(cos(t));for i=1:(length(t)-1)k1=F(t(i),y(i));k2=F(t(i)+h/2,y(i)+k1*h/2);k3=F(t(i)+h/2,y(i)+k2*h/2);k4=F(t(i)+h,y(i)+k3*h);y(i+1)=y(i)+1/6*(k1+2*k2+2*k3+k4)*h;endsubplot(211)plot(t,y,'-')xlabel('t');ylabel('y');title('Approximation');span=[0,10];p=y(1);[t1,y1]=ode45(F,span,p);subplot(212)plot(t1,y1)xlabel('t');ylabel('y');title('Exact');图1②一阶微分方程：()22*/x t x x t =- ，初始值x(1)=2，求解区间[1 3]。
MATLAB 程序： %%%%%%%%%%% 四阶龙哥库塔法解微分方程%%%%%%%%%%% x'(t)=(t*x-x^2)/t^2%%%%%%%%%%% x(1)=2, 1≤t ≤3, h=1/128%%%%%%%%%%% 精确解：x(t)=t/(0.5+lnt)h=1/128; %%%%% 步长tf=3;t=1:h:tf;x=zeros(1,length(t));x(1)=2; %%%%% 初始值F_tx=@(t,x)(t.*x-x.^2)./t.^2;for i=1:(length(t)-1)k_1=F_tx(t(i),x(i));k_2=F_tx(t(i)+0.5*h,x(i)+0.5*h*k_1);k_3=F_tx((t(i)+0.5*h),(x(i)+0.5*h*k_2));k_4=F_tx((t(i)+h),(x(i)+k_3*h));x(i+1)=x(i)+(1/6)*(k_1+2*k_2+2*k_3+k_4)*h; endsubplot(211)plot(t,x,'-');xlabel('t');ylabel('x');legend('Approximation');%%%%%%%%%%%%%%%%%%%%%%%%%%%% ode45求精确解t0=t(1);x0=x(1);xspan=[t0 tf];[x_ode45,y_ode45]=ode45(F_tx,xspan,x0);subplot(212)plot(x_ode45,y_ode45,'--');xlabel('t');ylabel('x');legend('Exact');图2二、四阶龙格库塔法解二阶微分方程①二阶微分方程：cos y t ，初始值y(0)=0，y'(0)=-1，求解区间[0 10]。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
MATLAB 在数值分析中的应用摘 要：数值分析是研究数学问题数值解及其理论的一个数学分支,涉及面很广.自计算机在其中应用以后,使其应用价值更为广泛,解决实际问题更为有效.以MA TLAB 为平台的数值分析方法与图形可视化编程在解决一些复杂问题时,大大降低了难度,减少了工作量.本文以经典龙格-库塔算法为核心,介绍用数值方法解常微分方程(组),分析自制系统的稳定性问题.如范德波尔方程.以及混沌现象中奇怪吸引子的模拟实现关键词： MA TLAB; 经典龙格-库塔算法; 解常微分方程(组); 范德波尔方程; 奇怪吸引子;1 MATLAB 平台1.1 矩阵MATLAB 是矩阵实验室（Matrix Laboratory ）的简称, 美国MathWorks 公司出品的商业数学软件.以下介绍在处理方程组时的特殊技法: A=[a b c; d e f; g h i]得到矩阵a b c d e f g h i也可直接输入方程组:function y=equition(t,x)y=[q(2)+q(1)/sqrt(q(1)^2+q(2)^2)*(1-q(1)^2-q(2)^2); % 以q 向量保存方程组;-q(1)+q(2)/sqrt(q(1)^2+q(2)^2)*(1-q(1)^2-q(2)^2)] % q(1)代表应变量x, q(2) 代表应变量y得到1.2 语言MATLAB 编程语言与 C 语言 类似, 格式上略有不同. 里不做详述, 程序中会做必要的标注.但在处理矩阵式十分方便, C++是以数组的形式存放变量,对某个元素都要准确指出其位置.22222222()]()]x y x y x y y x x y x y ⎧=+-+⎪+⎪⎨⎪=-+-+⎪+⎩C++是以数组的形式存放变量,对某个元素都要准确指出其位置,处理循环时算法复杂,即使使用指针也很容易出错.2经典龙格-库塔算法公式(15张)龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。
由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。
该算法是构建在数学支持的基础之上的。
对于一阶精度的欧拉公式有：yi＋１＝yi＋h*K1K1=f(xi,yi)当用点xi处的斜率近似值K1与右端点xi＋１处的斜率K2的算术平均值作为平均斜率K*的近似值，那么就会得到二阶精度的改进欧拉公式：yi＋１＝yi＋h*( K1+ K2)/2K1=f(xi,yi)K2=f(xi+h,yi+h*K1)依次类推，如果在区间[xi,xi+1]内多预估几个点上的斜率值K1、K2、……Km，并用他们的加权平均数作为平均斜率K*的近似值，显然能构造出具有很高精度的高阶计算公式。
经数学推导、求解，可以得出四阶龙格－库塔公式，也就是在工程中应用广泛的经典龙格－库塔算法：yi＋１＝yi＋h*( K1+ 2*K2 +2*K3+ K4)/6K1=f(xi,yi)K2=f(xi+h/2,yi+h*K1/2)K3=f(xi+h/2,yi+h*K2/2)K4=f(xi+h,yi+h*K3)3MATLAB 算法的实现3.1程序代码:function [tp yp] = RK45_n_h(dydt,tspan,y0,h)format long;tf = tspan(2);t = (tspan(1):h:tspan(2))'; % Bn=length(t); % Awhile t(n) < tft(n+1)=t(n)+hh; % 因为在t = (tspan(1):h:tspan(2))'累加过程中,n=n+1; % 可能出现大数吃小数的可能,这样可以确保积到上限if t(n) > tft(n)=tf; % A n依然指向积分上限break; % B t依然等间隔的存储着积分上下限endendtt=tspan(1);y(1,:)=y0;np=1; tp(np)=tt; yp(np,:)=y(1,:); %yp(np,:)用来记录函数值随时间的变化,i=1; %为了作图的精确方便,就以步长为间隔.while (1)tend = t(np + 1);hh = t(np + 1) - t(np);if hh > h, hh = h; endwhile(1)if tt + hh > tend, hh = tend - tt; endk1 = dydt(tt, y(i,:))';k2 = dydt(tt + hh/2, y(i,:) + k1.*hh./2)';k3 = dydt(tt + hh/2, y(i,:) + k2.*hh./2)';k4 = dydt(tt + hh, y(i,:) + k3*hh)';y(i+1,:) = y(i,:) + (k1 + 2*(k2 + k3) + k4)/6*hh;i = i + 1;tt=t(i);if tt >= tend, break, endendnp = np + 1; tp(np) = tt; yp(np,:) = y(i,:);if max(y(np,:))>1e200 %溢出时跳出循环break;endif tt >= tf, break, endendformat short3.2 应用力证function [tp yp] = RK45_n_h(dydt,tspan,y0,h)dydt是引用方程组, tspan 是积分上下限; y0是初值; h 为步长.函数调用格式:[A B]= RK45_n_h(@equition ,[0 ,10],[0.1, 0.1], 0.01)函数返回结果是两个向量tp 和 yp , tp 是从0开始,以0.01为步长,到10的向量; yp 是对应的函数值向量,有两个分量.初始值[0.1 0.1].输入命令:plot(B(:,1),B(:,2)),grid;得到实时解的平面相轨迹从图中反映出其解具有周期性,轨迹构成稳定的极限环.如果初值远离圆环 ,如[10 ,10], 解又是怎样的?[A B]= RK45_n_h(@equition ,[0 ,10],[10, 10], 0.01); plot(B(:,1),B(:,2));grid221x y +=得到解仍然是稳定的极限环为验证其准确性,可取时间跨度大一些[0 , 50],重复上述操作, 事实上还可以更直观的观察解的变化:plot(A,B(:,1),A,B(:,2)),grid可以看出随时间推移,x(1), x(2) 两个分量都是简谐振动, 相位相差 /2。
小结以上就是基于MA TLAB平台数值方法解常微分方程的大致过程4其他二维自治微分方程组的解法举例.4.1.1.a分别取不同的初值[1, 0.2]; [1, 0.5]; [10.8], [1，1.1]; [1, 1.2]axis([0 6 0 6]);xlabel('x轴');ylabel('y轴');[t,x]=ode23(ode913,[0 6],[1 0.2],option);[t,x]=ode23(ode913,[0 6],[1 0.5],option);[t,x]=ode23(ode913,[0 6],[1 0.8],option);[t,x]=ode23(ode913,[0 6],[1 1.1],option);[t,x]=ode23(ode913,[0 6],[1 1.2],option);hold onlegend('初值为[1,0.2]','初值为[1,0.5]','初值为[1,0.8]','初值为[1,1.1]','初值为[1,1.2]');text(4,0.5,'初值为[1,0.2]');text(3.5,2,'初值为[1,0.5]');text(2.2,2.5,'初值为[1,0.8]');text(1.3,3,'初值为[1,1.1]');text(0.9,2,'初值为[1,1.2]');得到 :不同初值,有不同的极限环系统是稳定的4.1.1.b 求如下系统的极限环，并判断其稳定性输入:[A B]= RK45_n_h(@BWD,[0 ,5],[0.001,0.00001], 0.001); plot(A,B(:,1),A,B(:,2)),grid 得到:[A B]= RK45_n_h(@BWD,[0 ,10],[0.01,0.1], 0.001); plot(A,B(:,1),A,B(:,2)),grid221/2222221/2222()(1)()(1)x x x y x y y y y x y x y x⎧=++-+⎪⎨=++--⎪⎩可见初值不同,只决定了开始是的轨迹,但在 的时候函数值也是趋于无穷的, (尽管在从 R=1 的圆内趋于圆的速度远小于从圆外趋于无穷的速度) 是不稳定的i.范德波尔方程 (解高阶常微分方程)将范德波尔方程化为状态方程当取 7时:编写M.函数function dy=verderpol (t,y) %定义输入,输出变量和函数文件名 dy-[y(2) ;7*(1-y(1)^2)*y(2)-y(1)]t →+∞2(1)0x x x x ε+-+=2(1)x yy x y x ε=⎧⎨=--⎩4.1.3 吸引子吸引子 能量耗散系统最终收缩到的一种定常状态。
这是一个动力系统在t →∞时所呈现的与时间无关的定态，并且不管选取什么样的初始值其终值的定态只有一个，也就是说终值与初始值无关。
这类吸引子也称平庸吸引子。
如：阻尼单摆有不动点吸引子，范德玻耳方程有极限环吸引子，等等。
奇怪吸引子 相对于平庸吸引子而言，它们的特点之一是终态值与初始值密切相关，或者说对初始值具有极端敏感性；初始取值的细微差别可能会导致完全不同的结果，这时的吸引子毫无周期可言，即所谓混沌。
4.1.3.a Lorenz 模型状态方程⎪⎪⎪⎪⎪⎨⎧--=--=dzxz y rx d dyy x d dxτστ)(在洛伦兹方程中，取参数s =10，b = 8/3，随参数r 增加，出现一次新分岔－霍夫分岔取r = 28 时计算的结果如下。
创建M.函数.function xdot=lorenzq(t,x)xdot=[-8/3*x(1)+x(2)*x(3);-10*x(2)+10*x(3);-x(1)*x(2)+28*x(2)-x(3)];在命令窗口输入[A B]= RK45_n_h(@lorenzq,[0 ,50],[0,0,1e-10], 0.001);plot(A,B),grid将 X 三个分量反映到的三维空间中得到下图plot3(B(:,1),B(:,2),B(:,3));grid不同时刻观察t=5t=50t=60从不同角度观察x-y平面x-z平面y-z平面奇怪吸引子的最重要特征是对初值的敏感性，初始相互靠近的两条轨线将按指数式规律分离。