为什么数值仿真里要用RK4(龙格库塔法)
小跳最近在搭建一个数值仿真环境,由于需要用到python里面的一些库,所以不得不把simulink的模型搬过来,我们都知道在simulink里,仿真的时候设置仿真步长和微分方程求解器是必要的步骤。但是为什么要设置这个小跳却早已忘记了。
一年级的时候搬砖搬多了,数分课也没好好上,回头一看,这么简单的东西,当时竟然整的稀里糊涂的。
为什么要用RK4
先po一张图,直观感受一下仿真的误差。
对于给定线性常微分方程 [dot x = x] 易得,其解是 [x(t) = Ce^t ] RK4是龙格库塔法曲线,None是一阶解法(x(t+dt) = x(t)+dot x dt) 可以看到,线性常微分方程误差尚且如此之大,那么推广到非线性微分方程,像这种形式 [ dot x = f(x,t) = tx^2 - frac{x}{t}... ] 那肯定误差直接起飞了。解析解求起来也挺麻烦,这里就不再引入分析了。
接下来把定义回顾一下,贴一下代码,有需自取,希望对大家有所帮助。
定义回顾
数值分析中,龙格-库塔法(Runge-Kutta methods)是用于非线性常微分方程的解的重要的一类隐式或显式迭代法。这些技术由数学家卡尔·龙格和马丁·威尔海姆·库塔于1900年左右发明。该方法主要是在已知方程导数和初值信息,利用计算机仿真时应用,省去求解微分方程的复杂过程。
令初值问题表述如下。 [ y' = f(t,y), y(t_0) = y_0 ] 则,对于该问题的RK4由如下方程给出: [ y_{n+1}=y_{n}+frac{h}{6}left(k_{1}+2 k_{2}+2 k_{3}+k_{4}right) \] 其中 [ begin{matrix} k_{1}=fleft(t_{n}, y_{n}right) \ k_{2}=fleft(t_{n}+frac{h}{2}, y_{n}+frac{h}{2} k_{1}right) \ k_{3}=fleft(t_{n}+frac{h}{2}, y_{n}+frac{h}{2} k_{2}right) \ k_{4}=fleft(t_{n}+h, y_{n}+h k_{3}right)end{matrix} ]
式中,(h)为仿真步长,满足(h<epsilon_1 rightarrow error<epsilon_2)
代码实现
import numpy as np
import matplotlib.pyplot as plt
import sympy as sy
# 这里介绍一个符号运算的方法,可以用来求解方程什么的
def diff_eq(t,x):
return sy.diff(x(t),t,1) - x(t)
t = sy.symbols('t')
x = sy.Function('x')
sy.pprint(sy.dsolve(diff_eq(t,x),x(t)))
def dot_x(t,x):
return x
def rk4(f,t,x,h):
k1 = f(t,x);
k2 = f(t+0.5*h,x + 0.5*h*k1)
k3 = f(t+0.5*h,x + 0.5*h*k2)
k4 = f(t+h,x + h*k3)
return h/6*(k1+2*k2+2*k3+k4)
t_list = np.arange(0,5,0.1);
#print(t)
x1_list = np.exp(t_list)
x2_list = []
x3_list = []
h = 0.1
x2 = 1;
x3 = 1;
for t in t_list:
# print(t,idx)
x2_list.append(x2)
x3_list.append(x3)
x2 = x2 + rk4(dot_x,t, x2, h)
x3 = x3 + dot_x(t,x3) * h
error_2 = x1_list - x2_list
error_3 = x1_list - x3_list
plt.figure()
plt.subplot(2,1,1)
plt.plot(t_list,x1_list, 'b-',label='Real')
plt.plot(t_list,x2_list,'r--', label = 'RK4')
plt.plot(t_list,x3_list,'g--', label = 'None')
plt.legend()
plt.subplot(2,1,2)
plt.plot(t_list,error_2, 'r--',label='Error_RK4')
plt.plot(t_list,error_3, 'g--',label='Error_none')
plt.legend()
plt.xlabel('Time(s)')
plt.show()
闲话
这里推荐一个提高效率的工具Matplotlib cheat sheet
对于一个经常画图的科研狗来说,这张图真是太太太太有必要了,因为时常遇到以下场景,不记得colormap名字,打开文档查一番,不记得线宽关键词,打开文档查一番,不记得marker名字,打开文档查一番。。。。。等等等等
所以,有了这张图,在平常画图的时候中遇到的95%需要查文档的问题都可以在这张图中找到答案。
这个速查表,可以关注微信公众号“探物及理”后台回复“python画图”领取。
- 深入源码探索 ReactNative 通信机制
- PHP跨站脚本攻击(XSS)漏洞修复思路(二)
- WordPress发布文章自动同步到新浪微博(带特色图片)
- go http 服务器编程(1)
- Linux系统内存监控、性能诊断工具vmstat命令详解
- go http 服务器编程(2)
- 利用placeholder属性来添加输入框默认文字提示,提高用户体验
- Linux系统监控、诊断工具之top命令详解
- 【Dev Club分享】iOS黑客技术大揭秘
- Linux终端:用cat命令查看不可见字符
- golang 函数定义及其接口实例
- 分享两种圣诞节雪花特效JS代码(网站下雪效果)
- React 移动 web 极致优化
- golang 高效字符串拼接
- JavaScript 教程
- JavaScript 编辑工具
- JavaScript 与HTML
- JavaScript 与Java
- JavaScript 数据结构
- JavaScript 基本数据类型
- JavaScript 特殊数据类型
- JavaScript 运算符
- JavaScript typeof 运算符
- JavaScript 表达式
- JavaScript 类型转换
- JavaScript 基本语法
- JavaScript 注释
- Javascript 基本处理流程
- Javascript 选择结构
- Javascript if 语句
- Javascript if 语句的嵌套
- Javascript switch 语句
- Javascript 循环结构
- Javascript 循环结构实例
- Javascript 跳转语句
- Javascript 控制语句总结
- Javascript 函数介绍
- Javascript 函数的定义
- Javascript 函数调用
- Javascript 几种特殊的函数
- JavaScript 内置函数简介
- Javascript eval() 函数
- Javascript isFinite() 函数
- Javascript isNaN() 函数
- parseInt() 与 parseFloat()
- escape() 与 unescape()
- Javascript 字符串介绍
- Javascript length属性
- javascript 字符串函数
- Javascript 日期对象简介
- Javascript 日期对象用途
- Date 对象属性和方法
- Javascript 数组是什么
- Javascript 创建数组
- Javascript 数组赋值与取值
- Javascript 数组属性和方法
- Thinkphp5框架ajax接口实现方法分析
- android MediaRecorder实现录屏时带录音功能
- php根据地址获取百度地图经纬度的实例方法
- Android 代码一键实现银行卡绑定功能
- Android 通过cmake的方式接入opencv的方法步骤
- Yii框架响应组件用法实例分析
- Android开发学习实现简单计算器
- Android Studio finish()方法的使用与解决app点击“返回”(直接退出)
- Android 8.1隐藏状态栏图标的实例代码
- Android制作登录页面并且记住账号密码功能的实现代码
- Yii框架分页技术实例分析
- PHP命名空间与自动加载机制的基础介绍
- Flutter下Android Studio配置gradle的方法
- Flutter 实现整个App变为灰色的方法示例
- Android studio开发小型对话机器人app(实例代码)