龙格库塔(runge-kutta，RK)法求解微分方程

求解微分方程的意思就是，已知导函数，求原函数。

先声明一点，欧拉法、中值法、龙哥库塔法求解微分方程，得出的结果不是表达式，而是一系列离散点。

一、欧拉法递推

问题描述：已知y'=f(x,y)，求y(x)。

例题：已知y'=y，y(0)=1，求y(x)。

解：这个题目高数知识很容易求解，答案是y=e^x。但是这里我们不使用解析法，而使用迭代递推法去求y(x)。

把y(x)在x0处泰勒展开到一阶，有y(x0+h)=y(x0)+y'(x0)h，有了这个式子，我们就能根据y(0)递推出y(0+h)的值，进而推出y(0+h+h)、y(0+h+h+h)的值····。根据这一原理就能很容易的写出递推程序了：

clear;
euler(0.1);%减小步长可以提高运算精度function euler(h)
%欧拉法解微分方程:y'=y     真实解为y=e^t
%形参：自变量的步长h
hold on;
idx = 1: 1: 10/h;x(1)=0;
y(1)=1;for i = idxx(i+1) = x(i) + h;dy = y(i);    y(i+1) = y(i) + h * dy;  %y1 = y0 + h * y0'
endplot(x, exp(x), 'r');
plot(x, y, 'b');
legend('真实解','递推解');
end

运行结果如下：

我们可以发现，递推的精度有点差，不过下面我们还有更好的递推方法。

二、中值法递推

该方法是对欧拉法的改进，原理就是拉格朗日中值定理，y(x0+h)=y(x0)+h*y'(a)，其中a∈[x0, x0+h]，虽然a不是区间的中点，但是这里我们认为用[x0, x0+h]区间中点处x0+h/2的导数来递推，比单纯用左端点的导数做递推，精度更高。该方法所用的递推公式为：

中值法有两种，第二种是这样的：用区间左右端点处的导数的平均值来递推，本质上就是左右端点各预测一个增量△y，这两个增量权重各占50%。

该方法的递推公式:

在接下来的章节中，我们将会看到，这两种中值方法，本质上都是二阶龙格库塔方法的特例。

下面是第一种中值法的程序：

clear;center(0.1);function center(h)
%中点法解微分方程:y'=y     真实解为y=e^t
%形参：自变量的步长h
hold on;idx = 1: 1: 10/h;x(1)=0;
y(1)=1;for i = idx    x(i+1) = x(i) + h;dy = y(i);    y_center = y(i) + h/2 * dy;%预测y(xi+h/2)函数值    dy = y_center;%中点的导数代替左端点的导数进行迭代y(i+1) = y(i) + h * dy;
endplot(x, exp(x), 'r');
plot(x, y, 'b');
legend('真实解','递推解');
end

可以发现，几乎重合，这个方法跟前面的欧拉法相比，步长相等的情况下，递推精度大幅提高，红蓝曲线几乎重合！

三、龙格库塔法（RungeKutta）求解

上面的方法2，看起来几乎重合，实际上放大以后还是能看到误差的，我们还有精度更高的方法。

https://wenku.baidu.com/view/d69e8f1f77c66137ee06eff9aef8941ea76e4b2c.html

推导过程会用到一元函数泰勒公式、复合函数求导公式。

泰勒公式：

复合函数求导公式（全导数公式）：

我们如果用泰勒二阶截断式来做递推，效果肯定比欧拉法要好，因为欧拉法的本质就是用的泰勒一阶截断式。

还是用前面的例题：已知y'=f(x,y)，y(0)=1，求y(x)。

用泰勒二阶递推公式能提高精度是显而易见的，但是由泰勒公式可见，其缺点是，通过提高阶次来运算时，需要用到二阶导、三阶导.....，有些高阶导并不是那么好求的，龙格库塔递推法就是对这一问题进行优化，用低阶导的泰勒展开来避免掉求高阶导。

龙格库塔(runge-kutta，RK)法求解微分方程相关推荐

龙格-库塔（Runge-Kutta）法解微分方程
龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法.由于此算法精度高,采取措施对误差进行抑制,所以其实现原理也较复杂.该算法是构建在数学支持的基础之上的. 对于一阶精度的欧 ...
阿当姆斯方法求解微分方程初值问题：四阶龙格库塔提供出发值（C语言）
#include<stdio.h>float f1(float x,float y) //第一问 {return y*y;//方程 }float f2(float x,float y)// ...
MATLAB从入门到精通-欧拉法与梯形法求解微分方程（含MATLAB源码）
前言以下是我为大家准备的几个精品专栏,喜欢的小伙伴可自行订阅,你的支持就是我不断更新的动力哟! MATLAB-30天带你从入门到精通 MATLAB深入理解高级教程(附源码) tableau可视化数据 ...
四阶龙格库塔法的基本思想_四阶龙格库塔实验报告.docx
四阶龙格库塔实验报告三.四阶Runge-Kutta法求解常微分方程一.龙格库塔法的思想根据第九章的知识可知道,Euler方法的局部截断误差是,而当用Euler方法估计出再用梯形公式进行校正,即采用改 ...
matlab求微分方程同届,Matlab学习——求解微分方程(组)
介绍: 1.在 Matlab 中,用大写字母 D 表示导数,Dy 表示 y 关于自变量的一阶导数,D2y 表示 y 关于自变量的二阶导数,依此类推.函数 dsolve 用来解决常微分方程(组)的求解问 ...
四阶龙格库塔（Runge-Kutta）求解微分方程-多种编程语言
前期是分享了matlab下面实现四阶龙格库塔(Runge-Kutta)求解微分方程,这期分享一下C++.C.Java.Python下面的四阶龙格库塔(Runge-Kutta)求解微分方程. 前文传送门 ...
matlab 龙格-库塔法求解常微分方程
最近学习分室模型,里面碰到了用matlab 龙格-库塔法求解常微分方程研究了一阵子终于明白到底怎么实现了: 1. matlab 新建.m文件,编写龙格-库塔法求解函数 function [x,y] ...
C#，数值计算，解微分方程的龙格-库塔二阶方法与源代码
微分方程含有导数或微分的方程称为微分方程,未知函数为一元函数的微分方程称为常微分方程. 微分方程的阶数微分方程中导数或微分的最高阶数称为微分方程的阶数. 微分方程的解使得微分方程成立的函数称为微 ...
[常微分方程的数值解法系列五] 龙格-库塔(RK4)法
龙格-库塔法简介基本思想具体方法一阶二阶求解参数特殊二阶三阶高阶步长选择例子在惯性导航以及VIO等实际问题中利用IMU求解位姿需要对IMU测量值进行积分得到需要的位置和姿态,其 ...

龙格库塔(runge-kutta，RK)法求解微分方程

龙格库塔(runge-kutta，RK)法求解微分方程相关推荐

最新文章

热门文章