四阶龙格库塔法求解微分方程【MATLAB】

2024-05-18 12:12:44

四阶龙格库塔法求解微分方程

求解过程数学描述

四阶龙格库塔的求解过程可用如下数学公式描述：
k1=f(tn,yn)k_1=f\left( t_n,y_n \right)k1=f(tn,yn)

k2=f(tn+h2,yn+h2k1)k_2=f\left( t_n+\frac{h}{2},y_n+\frac{h}{2}k_1 \right)k2=f(tn+2h,yn+2hk1)

k3=f(tn+h2,yn+h2k2)k_3=f\left( t_n+\frac{h}{2},y_n+\frac{h}{2}k_2 \right)k3=f(tn+2h,yn+2hk2)

k4=f(tn+h,yn+hk3)k_4=f\left( t_n+h,y_n+hk_3 \right)k4=f(tn+h,yn+hk3)

yn+1=yn+h6(k1+2k2+2k3+k4)y_{n+1}=y_n+\frac{h}{6}\left( k_1+2k_2+2k_3+k_4 \right)yn+1=yn+6h(k1+2k2+2k3+k4)

例子

为验证程序的有效性，选取一个已知解析解的微分方程验证。

y′=yy^{'}=yy′=y ，零初值状态，即y(0)=1y(0)=1y(0)=1，y=exy=e^xy=ex

code

编写的MATLAB程序如下：

% function RK4()
clc;clear;
Ts = 0.01;
h = Ts;
time = 1.5;
N = time/Ts;
t = linspace(Ts,time,N);y = zeros(1,N+1);
for m=2:Nk1 = exp(m*Ts);k2 = exp(m*Ts+h/2*k1);k3 = exp(m*Ts+h/2*k2);k4 = exp(m*Ts+h*k3);y(1,m+1) = y(1,m) +(k1+2*k2+2*k3+k4)*h/6;
end
figure
plot(t,exp(t))
hold on
y = y(1,2:N+1);
y = y+1;
plot(t,y)
legend('解析解','数值解');

结果分析

根据预期，算法应当逐渐收敛至稳态，但是实际的求解过程无法反应此过程。当求解区间变大时，出现算法出现轻微的发散现象，说明算法设计存在缺陷，原因尚未分析出来，后续理清思路后再补充。

参考文献

https://www.jianshu.com/p/e4aa9a688959

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

2022-10-26-四阶龙格库塔法计算程序修正

由于后续的工作需要使用数值计算方法，重写了四阶龙格库塔法，通过控制求解步长，可以有效的控制误差，上次遗留的发散问题仍未得到解决。

再次阅读之前写的程序，发现公布的算法是在反向验证龙格库塔算法的有效性，在解析解未知的前提下，算法无法进行正向求解。最新改进的算法已实现了正向求解。

%eqution
%y'=y y(0)=1clc;clear;% set the solution range and solution step
h = 0.01;
time = 5;
N = time/h;
t = linspace(h,time,N);
y = zeros(1,N);
y(1,1) = 1;% RK4
for m=1:N-1k1 = h*y(1,m);k2 = h*(y(1,m)+k1/2);k3 = h*(y(1,m)+k2/2);k4 = h*(y(1,m)+k3);y(1,m+1) = y(1,m) +(k1+2*k2+2*k3+k4)/6;
end% data visualization
figure
subplot(2,1,1);
plot(t,exp(t))
hold on
plot(t,y)
legend('解析解','数值解');
title('h=0.01')
subplot(2,1,2);
error = exp(t)-y;
plot(t,error)
legend('误差');

通过对比上方的两张图片，可以发现，在求解步长设置为0.01时，求解误差已经非常小。

四阶龙格库塔法求解微分方程【MATLAB】相关推荐

欧拉法、预估校正法(改进的欧拉法)与四阶龙格库塔法求解常微分方程的数值解C++程序
以y'=x+y,0<x<1,y(0)=1为例,取步长h=0.1,已知精确值为y=-x-1+2e^x,用来进行精度比较 #include<stdio.h> using names ...
四阶龙格库塔法求解一次常微分方程组(python实现)
四阶龙格库塔法求解一次常微分方程组一.前言二.RK4求解方程组的要点 1. 将方程组转化为RK4求解要求的标准形式 2. 注意区分每个方程的独立性三.python实现RK4求解一次常微分方程组 ...
龙格库塔法求解微分方程
在https://blog.csdn.net/weixin_42141390/article/details/110184743一文中,我们曾经讨论了欧拉法,龙格-库塔法也跟欧拉法一样,是用梯形的面积 ...
【Runge-Kutta】龙格-库塔法求解微分方程matlab仿真
1.软件版本 MATLAB2013b 2.算法理论龙格-库塔法(Runge-Kutta)是用于模拟常微分方程的解的重要的一类隐式或显式迭代法.龙格库塔法的家族中的一个成员如此常用,以至于经常被称为& ...
龙格库塔法和欧拉法求解微分方程的比较
文章目录计算机如何理解连续系统的动态特性? 欧拉法求解微分方程龙格库塔法求解微分方程 MATLAB代码编写和仿真效果计算机如何理解连续系统的动态特性? 一般连续系统的动态特性可以由一个微分方程, ...
用四阶龙格库塔法(RK4)求解二阶微分方程
import math import matplotlib.pyplot as plt# 龙格-库塔法的定义 def runge_kutta(y, h, f):k1 = h * f(t, x1, x2 ...
matlab解二阶微分方程组,[微分方程组]急急急！用MATLAB按二阶龙格库塔法求解微分方程组，急用于毕业设计！...
急急急!用MATLAB按二阶龙格库塔法求解微分方程组,急用于毕业设计! 问题补充:今天才发现自己之前做的一点都不对,17号就交论文了,我傻了,急死了!求各位大侠帮帮忙.谢谢!要求解的微分方程如图所示. ...
matlab龙格库塔法求通解,基于matlab及龙格库塔法求解布拉修斯方程.doc
基于matlab及龙格库塔法求解布拉修斯方程 Runge-Kutta法求解布拉修斯解摘要薄剪切层方程主要有三种解法,即相似解,非相似条件下对偏微分方程组的数值解和近似解.布拉修斯解是布拉修斯于19 ...
四阶龙格库塔法的基本思想_利用龙格库塔法求解郎之万方程.doc
利用龙格库塔法求解郎之万方程.doc 利用龙格-库塔法求解朗之万方程1. 待解问题布朗颗粒是非常微小的宏观颗粒,其直径的典型大小为10-710-6m.颗粒不断受到液体介质分子的碰撞,在任一瞬间,一个颗 ...
Matlab求解微分方程数值解
有三种方法求解微分方程数值解: 欧拉法改进欧拉法龙格库塔法接下来用一个练习来对比这三种求解方法. 问题描述:用改进的Euler方法.MATLAB的ode45命令分别求下列初值问题的数值解,并画图 ...

最新文章

热门文章