龙格-库塔法求微分方程2

Python是一种高级编程语言，广泛应用于科学计算和工程领域。

而在科学计算中，求解二阶常微分方程是一个常见的问题。

龙格库塔法（Runge-Kutta method）是一种常用的数值求解方法，可以用来求解常微分方程的初值问题。

在Python中，我们可以利用其强大的数值计算库来实现龙格库塔法求解二阶常微分方程。

在本文中，我们将介绍如何使用Python中的龙格库塔法来求解二阶常微分方程。

文章将从以下几个方面展开讲解：1. 二阶常微分方程的基本概念2. 龙格库塔法的原理与公式推导3. Python中如何实现龙格库塔法求解二阶常微分方程4. 一个具体的求解例子及代码实现5. 总结与展望一、二阶常微分方程的基本概念二阶常微分方程是指具有形如y''(t) = f(t, y, y')(t)的形式，其中t是自变量，y是未知函数，f是已知函数。

求解二阶常微分方程的目标是找到一个满足方程的函数y(t)。

通常情况下，需要给出该方程的初值条件，即y(0)和y'(0)，以求得方程的具体解。

二、龙格库塔法的原理与公式推导龙格库塔法是一种数值求解常微分方程初值问题的方法，通过迭代计算来逼近方程的解。

它的基本思想是将微分方程的解按照一定的步长进行逼近，然后利用逼近值来计算下一个逼近值，直到达到所需的精度。

龙格库塔法的一般形式为：k1 = h * f(tn, yn)k2 = h * f(tn + 1/2 * h, yn + 1/2 * k1)k3 = h * f(tn + 1/2 * h, yn + 1/2 * k2)k4 = h * f(tn + h, yn + k3)yn+1 = yn + 1/6 * (k1 + 2*k2 + 2*k3 + k4)其中，h是步长，t是自变量，y是未知函数，f是已知函数，k1、k2、k3、k4是中间变量。

三、Python中如何实现龙格库塔法求解二阶常微分方程在Python中，可以利用其强大的数值计算库来实现龙格库塔法求解二阶常微分方程。

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函数求解该一阶常微分方程组即可得到原二元二阶常微分方程组的数值解。

龙格库塔求解二阶常微分方程一、前言在数学领域中，常微分方程是一种非常重要的数学工具。

它可以用来描述许多自然现象，如物理学、化学、生物学等。

而龙格库塔法是一种求解常微分方程的方法之一。

本文将介绍如何使用龙格库塔法求解二阶常微分方程。

二、什么是二阶常微分方程二阶常微分方程是指形如y''+p(t)y'+q(t)y=f(t)的方程，其中y表示未知函数，p(t)、q(t)和f(t)都是已知函数。

通常情况下，我们需要找到一个特定的y函数来满足这个方程。

三、什么是龙格库塔法龙格库塔法是一种数值求解常微分方程的方法。

它基于欧拉方法，但更准确和稳定。

欧拉方法使用线性逼近来估计未知函数值，而龙格库塔法使用高阶多项式逼近来更准确地估计未知函数值。

四、龙格库塔法的基本思想1. 将时间区间[0,T]平均分成N个小区间；2. 在每个小区间内进行递推计算，直到得到所有时间点上的y值；3. 每次递推计算都使用龙格库塔公式来估算y的值。

五、龙格库塔法的公式对于二阶常微分方程y''+p(t)y'+q(t)y=f(t)，我们可以使用以下公式来递推计算：1. k1=h*y'(t)l1=h*f(t)-p(t)*k1/2-q(t)*y(t);2. k2=h*(y'(t)+l1/2)l2=h*f(t+h/2)-p(t+h/2)*k2/2-q(t+h/2)*(y(t)+k1/2);3. k3=h*(y'(t)+l2/2)l3=h*f(t+h/2)-p(t+h/2)*k3/2-q(t+h/2)*(y(t)+k2/2);4. k4=h*(y'(t)+l3)l4=h*f(t+h)-p(t+h)*k4-q(t+h)*(y(t)+k3);其中，h表示时间步长，t表示当前时间点，f表示已知函数，p和q都是已知常数。

通过这些公式，我们可以得到下一个时间点上的y值。

六、龙格库塔法的实现为了更好地理解龙格库塔法，我们可以看一下具体的实现过程。

二阶微分方程龙格库塔法
1.什么是二阶微分方程？
二阶微分方程是指二阶导数出现的方程，其通用形式为
y''+P(x)y'+Q(x)y=F(x)，其中P(x)、Q(x)、F(x)均为已知函数，y是所求函数。

2.为什么要用龙格库塔法？
求解二阶微分方程是数学、物理等领域中常见的问题，然而大多数情况下无法直接解题，所以需要使用数值方法。

龙格库塔法是一种数值方法，被广泛应用于求解微分方程，其优点是精度高、计算速度快。

3.龙格库塔法的原理
龙格库塔法是一种迭代方法，将微分方程看作初值问题，从初始点出发，采用一定的步长h，在每个点上用插值公式逼近y(x+h)，以此求得y(x+h)的近似值，一步步逼近所要求的精度。

4.龙格库塔法的步骤
（1）确定步长h和积分区间[a,b]。

（2）用初值y(x0)=y0，y'(x0)=y'0求出y(x0+h)的近似值。

（3）根据龙格库塔公式求得y(x0+2h)的近似值。

（4）对于连续求解的情况，重复以上步骤，直到求得所需的精度或者达到指定的终点。

5.龙格库塔公式
龙格库塔法的精度与所采用的公式有关，一般采用二阶或四阶的龙格库塔公式。

二阶龙格库塔公式为：
y0=y(x0)
k1=h*f(x0,y0)
k2=h*f(x0+h,y0+k1)
y(x0+2h)=y0+1/2(k1+k2)
其中，f(x,y)是待求函数。

6.总结
龙格库塔法是求解微分方程的一种常用数值方法，可以高精度、高效地解决二阶微分方程的问题。

该方法所需的计算量较小，易于编写程序实现，在实际应用中具有广泛的用途。

四阶龙格库塔解二阶微分方程四阶龙格库塔法是一种常用的数值解微分方程的算法，适用于解一阶、二阶以及高阶微分方程。

本文将围绕“四阶龙格库塔解二阶微分方程”展开，并分步骤进行阐述。

第一步：转换为一阶微分方程组对于二阶微分方程 y''=f(x,y,y')，可以将其转换为一阶微分方程组。

令 u=y'，则 y''=du/dx，原方程变为 du/dx=f(x,y,u)，整理得到以下一阶微分方程组：dy/dx=udu/dx=f(x,y,u)第二步：确定步长和初始条件在使用四阶龙格库塔法时，需要确定一个步长 h，以及初始条件y0 和 u0。

步长 h 决定了每一步积分时横坐标 x 的变化量。

初始条件 y0 和 u0 分别代表了 x=0 时的纵坐标和导数值。

第三步：进行数值积分按照四阶龙格库塔法的步骤进行数值积分，可得到第 k 步的纵坐标 yk 和导数值 uk：k1 = hf(xk,yk,uk)j1 = hukk2 = hf(xk+h/2,yk+j1/2,uk+k1/2)j2 = h(uk+k1/2)k3 = hf(xk+h/2,yk+j2/2,uk+k2/2)j3 = h(uk+k2/2)k4 = hf(xk+h,yk+j3,uk+k3)j4 = h(uk+k3)yk+1 = yk+(j1+2j2+2j3+j4)/6uk+1 = uk+(k1+2k2+2k3+k4)/6其中，k1 到 k4 和 j1 到 j4 分别为由函数 f(x,y,u) 求得的中间变量。

yk 和 uk 即为求得的第 k 步的纵坐标和导数值，yk+1 和uk+1 分别为求得的第 k+1 步的纵坐标和导数值。

第四步：循环迭代将上述计算步骤循环执行至所需的精度或区间内积分完成为止。

即对于每一步 k，根据当前的 yk 和 uk 求解下一步的 yk+1 和 uk+1，并记录下这些变量的值。

至此，关于“四阶龙格库塔解二阶微分方程”的阐述就结束了。

四阶龙格库塔方法使用四级四阶经典显式Rungkutta公式稳定性很好，RK4法是四阶方法，每步的误差是h5阶，而总积累误差为h4阶。

所以比欧拉稳定。

运行第三个程序：在一幅图中显示欧拉法和RK4法，随着截断误差的积累，欧拉法产生了较大的误差h=0.01h=0.0001，仍然是开始较为稳定，逐渐误差变大总结：RK4是很好的方法，很稳定，而且四阶是很常用的方法，因为到五阶的时候精度并没有相应提升。

通过这两种方法计算出角度峰值y=3.141593，周期是1.777510。

三个程序欧拉法clear;clch=0.00001;a=0;b=25;x=a:h:b;y(1)=0;z(1)=0;for i=1:length(x)-1 % 欧拉y(i+1)=y(i)+h*z(i);z(i+1)=z(i)+h*7.35499*cos(y(i));endplot(x,y,'r*');xlabel('时间');ylabel('角度');A=[x,y];%y(find(y==max(y)))%Num=(find(y==max(y)))[y,T]=max(y);fprintf('角度峰值等于%d',y) %角度的峰值也就是πfprintf('\n')fprintf('周期等于%d',T*h)%周期legend('欧拉');龙格库塔方法先定义函数rightf_sys1.mfunction w=rightf_sys1(x,y,z)w=7.35499*cos(y);clear;clc;%set(0,'RecursionLimit',500)h=0.01;a=0;b=25;x=a:h:b;RK_y(1)=0; %初值%RK_z(1)=0;初值for i=1:length(x)-1K1=RK_z(i); L1=rightf_sys1(x(i),RK_y(i),RK_z(i));%K1 and L1K2=RK_z(i)+0.5*h*L1;L2=rightf_sys1(x(i)+0.5*h,RK_y(i)+0.5*h*K1,RK_z(i)+0.5*h*L1);K3=RK_z(i)+0.5*h*L2;L3=rightf_sys1(x(i)+0.5*h,RK_y(i)+0.5*h*K2,RK_z(i)+0.5*h*L2);K4=RK_z(i)+h*L3;% K4L4=rightf_sys1(x(i)+h,RK_y(i)+h*K3,RK_z(i)+h*L3);and L4RK_y(i+1)=RK_y(i)+1/6*h*(K1+2*K2+2*K3+K4);RK_z(i+1)=RK_z(i)+1/6*h*(L1+2*L2+2*L3+L4);endplot(x,RK_y,'b+');xlabel('Variable x');ylabel('Variable y');A=[x,RK_y];[y,T]=max(RK_y);legend('RK4方法');fprintf('角度峰值等于%d',y) %角度的峰值也就是πfprintf('\n')%周期fprintf('周期等于%d',T*h)两个方法在一起对比使用跟上一个相同的函数rightf_sys1.mclear;clc; %清屏h=0.0001;a=0;b=25;x=a:h:b;Euler_y(1)=0;%欧拉的初值Euler_z(1)=0;RK_y(1)=0;%龙格库塔初值RK_z(1)=0;for i=1:length(x)-1%先是欧拉法Euler_y(i+1)=Euler_y(i)+h*Euler_z(i);Euler_z(i+1)=Euler_z(i)+h*7.35499*cos(Euler_y(i));%龙格库塔K1=RK_z(i); L1=rightf_sys1(x(i),RK_y(i),RK_z(i)); % K1 andL1K2=RK_z(i)+0.5*h*L1;L2=rightf_sys1(x(i)+0.5*h,RK_y(i)+0.5*h*K1,RK_z(i)+0.5*h*L1);% K2 and L2K3=RK_z(i)+0.5*h*L2;L3=rightf_sys1(x(i)+0.5*h,RK_y(i)+0.5*h*K2,RK_z(i)+0.5*h*L2);% K3 and L3K4=RK_z(i)+h*L3; L4=rightf_sys1(x(i)+h,RK_y(i)+h*K3,RK_z(i)+h*L3); K4 and L4RK_y(i+1)=RK_y(i)+1/6*h*(K1+2*K2+2*K3+K4);RK_z(i+1)=RK_z(i)+1/6*h*(L1+2*L2+2*L3+L4);endplot(x,Euler_y,'r-',x,RK_y,'b-');[y,T]=max(RK_y);%角度的峰值也就是πfprintf('角度峰值等于%d',y)fprintf('\n')%周期fprintf('周期等于%d',T*h)xlabel('时间');ylabel('角度');legend('欧拉','RK4');。

四阶龙格库塔方法求解n自由度二阶微分方程在数值计算中，龙格-库塔方法（Runge-Kutta method）是一种常用的求解常微分方程（ODE）的数值方法。

它是一种多步法（multiple-step method），通过使用不同的近似值来迭代计算ODE 的解。

其中，四阶龙格-库塔方法是最常用的龙格-库塔方法之一，也是一种高阶方法，可用于求解n自由度二阶微分方程。

下面将介绍四阶龙格-库塔方法的原理、步骤和应用。

首先，我们来回顾一下二阶微分方程的一般形式：y''(y)=y(y,y(y),y'(y))这里，y(y,y(y),y'(y))是已知的函数，y(y)是我们要求解的未知函数。

为了求解这个二阶微分方程，我们需要将其转化为一个一阶微分方程的求解问题。

令y(y)=y(y)y(y)=y'(y)这样，我们可以将原始的二阶微分方程转化为如下的一阶微分方程组：y'(y)=y(y)y'(y)=y(y,y(y),y(y))现在，我们可以利用四阶龙格-库塔方法来求解这个一阶微分方程组。

四阶龙格-库塔方法基于泰勒展开的思想，通过使用一系列的近似值来逼近方程的解。

该方法的一般形式为：y1=ℎy(yy,yy)y2=ℎy(yy+ℎ/2,yy+y1/2)y3=ℎy(yy+ℎ/2,yy+y2/2)y4=ℎy(yy+ℎ,yy+y3)yy+1=yy+1/6(y1+2y2+2y3+y4)yy+1=yy+ℎ其中，yy是当前时间步，yy+1是下一个时间步，yy是当前位置的近似解，yy+1是下一个位置的近似解，ℎ是时间步长，y1、y2、y3和y4是中间的近似值。

四阶龙格-库塔方法的步骤如下：Step 1:给定初始条件我们需要提供初始条件，包括初始位置y0和初始速度y0，以及时间步长ℎ、求解的时间区间[y0,yy]和离散时间步数y。

y0=y(y0)y0=y'(y0)Step 2:进行迭代计算从n=0开始，进行y步的迭代计算，其中y=(yy−y0)/ℎ。

数值计算：四阶龙格-库塔法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%。

二阶龙格库塔法例题二阶龙格-库塔法（RK2）是一种数值求解常微分方程（ODE）的方法，它基于欧拉法的改进。

下面我将以一个例题来说明二阶龙格-库塔法的应用。

假设我们要求解如下的一阶常微分方程：dy/dx = x + y.初始条件为 y(0) = 1。

我们希望使用二阶龙格-库塔法来计算在 x = 0.1 处的函数值。

首先，我们需要将方程转化为标准形式，即将一阶微分方程表示为 dy/dx = f(x, y) 的形式。

在这个例子中，f(x, y) = x + y。

接下来，我们将计算步长 h，这是我们在每个步骤中前进的距离。

在二阶龙格-库塔法中，我们需要选择一个合适的步长。

通常情况下，我们可以通过试验不同的步长来选择一个合适的值。

在这个例子中，我们选择 h = 0.1。

然后，我们可以开始应用二阶龙格-库塔法的迭代公式。

迭代公式如下：k1 = h f(x, y)。

k2 = h f(x + h/2, y + k1/2)。

y_new = y + k2。

现在，我们可以开始迭代计算。

根据初始条件，我们有 x = 0，y = 1。

我们将进行以下计算步骤：1. 计算 k1 = h f(x, y) = 0.1 (0 + 1) = 0.1。

2. 计算 k2 = h f(x + h/2, y + k1/2) = 0.1 (0.05 + 1 + 0.1/2) = 0.105。

3. 计算 y_new = y + k2 = 1 + 0.105 = 1.105。

这样，我们得到了在 x = 0.1 处的函数值 y_new = 1.105。

通过不断重复上述步骤，我们可以计算出在其他点的函数值。

只需将当前的 x 和 y_new 作为下一步的初始条件，重复上述计算步骤即可。

需要注意的是，二阶龙格-库塔法是一种数值近似方法，其结果并不完全精确。

因此，在实际应用中，我们需要根据问题的要求和精度需求选择合适的步长和迭代次数。

以上就是使用二阶龙格-库塔法求解一阶常微分方程的例子。

隆格库塔法求解常微分方程摘要科学技术中常常需要求解常微分方程的定解问题，这里问题最简单的形式，是本章将着重考察的一阶方程的初值问题.虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些特殊类型的方程，实际问题中归结出来的微分方程主要靠数值解法求解.本文着重讨论了隆格库塔法求解一阶常微分方程的初值问题，采用了精度较高的经典的四阶隆格库塔法,然后通过对实例运用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)。

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。

Matlab 应用使用Euler 和Rungkutta 方法解臂状摆的能量方程背景 单摆是常见的物理模型，为了得到摆角θ的关于时间的函数，来描述单摆运动。

由角动量定理我们知道εJ M =化简得到0sin 22=+θθl g dt d 在一般的应用和计算中，只考虑摆角在5度以内的小摆动，因为可以吧简化为，这样比较容易解。

实际上这是一个解二阶常微分方程的问题。

θ在这里的单摆是一种特别的单摆，具有均匀的质量M 分布在长为2的臂状摆上，使用能量法建立方程W T =h mg ∆=2J 21ω化简得到四阶龙格库塔方法使用四级四阶经典显式Rungkutta公式稳定性很好，RK4法是四阶方法，每步的误差是h5阶，而总积累误差为h4阶。

所以比欧拉稳定。

运行第三个程序：在一幅图中显示欧拉法和RK4法，随着截断误差的积累，欧拉法产生了较大的误差h=0.01h=0.0001，仍然是开始较为稳定，逐渐误差变大总结：RK4是很好的方法，很稳定，而且四阶是很常用的方法，因为到五阶的时候精度并没有相应提升。

通过这两种方法计算出角度峰值y=3.141593，周期是1.777510。

三个程序欧拉法clear;clch=0.00001;a=0;b=25;x=a:h:b;y(1)=0;z(1)=0;for i=1:length(x)-1 % 欧拉y(i+1)=y(i)+h*z(i);z(i+1)=z(i)+h*7.35499*cos(y(i));endplot(x,y,'r*');xlabel('时间');ylabel('角度');A=[x,y];%y(find(y==max(y)))%Num=(find(y==max(y)))[y,T]=max(y);fprintf('角度峰值等于%d',y) %角度的峰值也就是πfprintf('\n')fprintf('周期等于%d',T*h)%周期legend('欧拉');龙格库塔方法先定义函数rightf_sys1.mfunction w=rightf_sys1(x,y,z)w=7.35499*cos(y);clear;clc;%set(0,'RecursionLimit',500)h=0.01;a=0;b=25;x=a:h:b;RK_y(1)=0; %初值%RK_z(1)=0;初值for i=1:length(x)-1K1=RK_z(i); L1=rightf_sys1(x(i),RK_y(i),RK_z(i));%K1 and L1K2=RK_z(i)+0.5*h*L1;L2=rightf_sys1(x(i)+0.5*h,RK_y(i)+0.5*h*K1,RK_z(i)+0.5*h*L1);K3=RK_z(i)+0.5*h*L2;L3=rightf_sys1(x(i)+0.5*h,RK_y(i)+0.5*h*K2,RK_z(i)+0.5*h*L2);K4=RK_z(i)+h*L3;% K4L4=rightf_sys1(x(i)+h,RK_y(i)+h*K3,RK_z(i)+h*L3);and L4RK_y(i+1)=RK_y(i)+1/6*h*(K1+2*K2+2*K3+K4);RK_z(i+1)=RK_z(i)+1/6*h*(L1+2*L2+2*L3+L4);endplot(x,RK_y,'b+');xlabel('Variable x');ylabel('Variable y');A=[x,RK_y];[y,T]=max(RK_y);legend('RK4方法');fprintf('角度峰值等于%d',y) %角度的峰值也就是πfprintf('\n')%周期fprintf('周期等于%d',T*h)两个方法在一起对比使用跟上一个相同的函数rightf_sys1.mclear;clc; %清屏h=0.0001;a=0;b=25;x=a:h:b;Euler_y(1)=0;%欧拉的初值Euler_z(1)=0;RK_y(1)=0;%龙格库塔初值RK_z(1)=0;for i=1:length(x)-1%先是欧拉法Euler_y(i+1)=Euler_y(i)+h*Euler_z(i);Euler_z(i+1)=Euler_z(i)+h*7.35499*cos(Euler_y(i));%龙格库塔K1=RK_z(i); L1=rightf_sys1(x(i),RK_y(i),RK_z(i)); % K1 andL1K2=RK_z(i)+0.5*h*L1;L2=rightf_sys1(x(i)+0.5*h,RK_y(i)+0.5*h*K1,RK_z(i)+0.5*h*L1);% K2 and L2K3=RK_z(i)+0.5*h*L2;L3=rightf_sys1(x(i)+0.5*h,RK_y(i)+0.5*h*K2,RK_z(i)+0.5*h*L2);% K3 and L3K4=RK_z(i)+h*L3; L4=rightf_sys1(x(i)+h,RK_y(i)+h*K3,RK_z(i)+h*L3); K4 and L4RK_y(i+1)=RK_y(i)+1/6*h*(K1+2*K2+2*K3+K4);RK_z(i+1)=RK_z(i)+1/6*h*(L1+2*L2+2*L3+L4);endplot(x,Euler_y,'r-',x,RK_y,'b-');[y,T]=max(RK_y);%角度的峰值也就是πfprintf('角度峰值等于%d',y)fprintf('\n')%周期fprintf('周期等于%d',T*h)xlabel('时间');ylabel('角度');legend('欧拉','RK4');。