四阶龙格库塔法的基本思想_数值常微分方程-欧拉法与龙格-库塔法
大三时候在跳蚤市场闲逛,从一位数学院的学长那里买了一些闲书,最近翻出来刚好有李荣华、刘播老师的《微分方程数值解法》和王仁宏老师的《数值逼近》,结合周善贵老师的《计算物理》课程,整理一下笔记。
本文整理常微分方程数值求解的欧拉法与龙格-库塔法。
一般地,动力学系统的时间演化可以用常微分方程的初值问题来描述,例如设一维简谐运动的回复力:
因此本文主要整理一阶常微分方程初值问题的数值解法。
一阶常微分方程初值问题
设
存在常数L,使得
假设
欧拉法
将区间
推导
1、根据泰勒展开式:
略去二阶小量,得:
以此类推,得到递推公式:
2、数值积分推导
由
以此类推,可得到:
为了提高精度,可以使用梯形积分代替矩形积分,即:
以此类推,得到改进的欧拉法:
Python计算实例
以
import math
from matplotlib import pyplot as pltt_0 = 0
y_0 = 1
tau = 0.1
i = 1
solve = []
Euler = []
t = []
while i < 100:if i == 1:y_n = y_0t_n = t_0Euler.append(y_n)solve.append(math.exp(t_n))t.append(t_n)func = y_ny_n = y_n + tau * funct_n = t_n + taui += 1plt.plot(t, Euler, c='green', label=' Euler method')
plt.plot(t, solve, c='red', label=' accuracy')
plt.fill_between(t, solve, Euler, facecolor='blue', alpha=0.2)
plt.title('Euler method', fontsize=19)
plt.xlabel('t', fontsize=19)
plt.ylabel('y', fontsize=19)
plt.legend()
plt.show()
作图可以看到,当迭代步数较多后,欧拉法的结果逐渐落后于精确指数解的增长速度。下面分析欧拉法的误差来源。
在
在计算中,我们更关心精确解和数值解之间的误差
根据Lipschitz条件,可得:
由
局部截断误差
稳定性分析
如果计算的初值不能精确给定,例如存在测量、舍入误差等,在计算过程中,每一步传递的误差连续依赖于初始误差,则称算法稳定,否则该算法不稳定。
对于不同的初值
两式相减,得:
根据Lipschitz条件,可得:
龙格-库塔法
龙格库塔法的主要思想:在
将
令
以此类推,可以得到:
同时,我们可以写出泰勒展开的形式解:
其中:
通项为:
基本思路是,利用当前点的函数值
现在把
将
将
与泰勒展开式
2个方程有3个未知数,因此有无穷多个解,可采用
令
此即为二阶龙格-库塔法。
与上一节的欧拉法公式对比:
Python计算实例
仍以
import math
from matplotlib import pyplot as pltt_0 = 0
y_0 = 1
z_0 = 1
tau = 0.1
i = 1
j = 1
solve = []
Euler = []
R_K = []
t = []
while i < 100:if i == 1:y_n = y_0t_n = t_0R_K.append(y_n)solve.append(math.exp(t_n))t.append(t_n)func_n = y_nfunc_m = y_n + tau * func_ny_n = y_n + 0.5 * tau * (func_n + func_m)t_n = t_n + taui += 1
t = []
while j < 100:if j == 1:z_n = z_0t_n = t_0Euler.append(z_n)t.append(t_n)func = z_nz_n = z_n + tau * funct_n = t_n + tauj += 1plt.scatter(t, R_K, marker='^', c='blue', s=70, label=' R-K method')
plt.plot(t, Euler, c='green', label=' Euler method')
plt.plot(t, solve, c='red', label=' accuracy')
plt.fill_between(t, solve, Euler, facecolor='yellow', alpha=0.2)
plt.title('Euler method & R-K method', fontsize=19)
plt.xlabel('t', fontsize=19)
plt.ylabel('y', fontsize=19)
plt.legend()
plt.show()
黄色部分表示数值解和精确解的偏离,可以看到,二阶龙格-库塔法(改进的欧拉法)精确度得到了很大的提升。
二阶龙格-库塔法中,泰勒展开到了
Reference:
1、周善贵,《计算物理》课程讲义
2、李荣华,刘播,《微分方程数值解法》
3、王仁宏,《数值逼近》
四阶龙格库塔法的基本思想_数值常微分方程-欧拉法与龙格-库塔法相关推荐
- python解常微分方程龙格库_数值常微分方程-欧拉法与龙格-库塔法
大三时候在跳蚤市场闲逛,从一位数学院的学长那里买了一些闲书,最近翻出来刚好有李荣华.刘播老师的<微分方程数值解法>和王仁宏老师的<数值逼近>,结合周善贵老师的<计算物理& ...
- 二阶龙格库塔公式推导_数值常微分方程-欧拉法与龙格-库塔法
大三时候在跳蚤市场闲逛,从一位数学院的学长那里买了一些闲书,最近翻出来刚好有李荣华.刘播老师的<微分方程数值解法>和王仁宏老师的<数值逼近>,结合周善贵老师的<计算物理& ...
- 四阶龙格库塔法的基本思想_利用龙格库塔法求解郎之万方程.doc
利用龙格库塔法求解郎之万方程.doc 利用龙格-库塔法求解朗之万方程1. 待解问题布朗颗粒是非常微小的宏观颗粒,其直径的典型大小为10-710-6m.颗粒不断受到液体介质分子的碰撞,在任一瞬间,一个颗 ...
- 四阶龙格库塔法的基本思想_四阶龙格库塔实验报告.docx
四阶龙格库塔实验报告 三.四阶Runge-Kutta法求解常微分方程一.龙格库塔法的思想根据第九章的知识可知道,Euler方法的局部截断误差是,而当用Euler方法估计出再用梯形公式进行校正,即采用改 ...
- 四阶龙格库塔法的基本思想_经典四阶龙格库塔法解一阶微分方程组讲义.doc
1.经典四阶龙格库塔法解一阶微分方程组 1.1运用四阶龙格库塔法解一阶微分方程组算法分析 , 经过循环计算由 推得 -- 每个龙格-库塔方法都是由一个合适的泰勒方法推导而来,使得其最终全局误差为,一种 ...
- 四阶龙格库塔法的基本思想_龙格库塔积分算法
龙格库塔法 龙格库塔法是常用于模拟常微分方程的解的重要的一类隐式或显式迭代法.这些技术由数学家C. Runge和M.W. Kutta于1900年左右发明.由于此算法精度高,采取措施对误差进行抑制,所以 ...
- 四阶龙格库塔法的基本思想_请问用四阶龙格库塔法解二阶微分方程的思想是什么?...
默认y的单位是弧度 k=1000; t=0:0.001:1; Y=[]; err=1 K=[]; Ymax=[]; xishu=1.01; while err X=[0 0]; k=xishu*k; ...
- 四阶龙格库塔法的基本思想_Runge-Kutta法求四元数微分方程
Runge-Kutta法求四元数微分方程 Runge-Kutta法求四元数微分方程 文章目录一.背景知识1. 坐标系 2. 四元数四元数的矩阵形式 四元数与旋转的关系 二.数学模型1. 四元数微分方程 ...
- 常微分方程组之龙格-库塔法
对于方程y'=f(x,y),初始条件:y(x0)=y0, #include<stdio.h>FILE *fp=fopen("ex的值.dat","w" ...
最新文章
- java运行按钮在哪里_[tkinter按钮命令已在程序启动时运行
- 测开之路二十:比较v1和v2
- IDEA和VS code设置默认换行符为LF
- 查看Cglib生成的Class(字节码)文件
- Linux下dig命令使用
- 综述: 通信雷达一体化中的信号处理
- ElementUI表格表头对角线的绘制
- Word中大括号内公式左对齐
- k8s deployment Strategy 更新策略
- Error “Client wants topic A to have B, but our version has C. Dropping connection.“
- 2019年暑假第八周总结
- html页面国际化之谷歌翻译js实践,支持通过判断浏览器语言自动将中文翻译成英文
- 陆源:阿贝尔对椭圆函数论的贡献[附椭圆函数、模形式(g_2,g_3)、模函数的C++程序计算]
- Github搜索开源项目过滤技巧
- swift class的虚函数表、扩展、@objc修饰、虚函数的派发方式研究
- const与指针用法
- 简单登录功能(一)token的使用
- php 导致服务器成肉鸡,把我的服务器搞成了ddos肉鸡,如何解决?
- ser2net的测试
- 计算机学院校运会解说词,运动会学院解说词
热门文章
- 织梦其他模型使用联动类型地区联动
- HDU - 1024 Max Sum Plus Plus 最大m段子段和+滚动数组优化
- 选哪个云计算平台部署自己的网站?
- CSS Hack 汇总快查
- 学习网页栅格系统的几篇好文
- 鸿蒙3部曲先看哪部,讨论雪鹰与鸿蒙三部曲的关系
- php curl 采集文件,curl获取远程文件内容
- es springboot 不设置id_springboot整合ES_文档ID删除
- 福州java培训哪里好_南通java培训哪家好
- android 网易item广告,Android仿网易严选商品详情页