常微分方程初值问题数值解法[完整公式]（Python）
2024-04-13 13:07:53
目录
1、概述
（1）常微分方程初值问题数值解法
（2）解题步骤
（3）数值微分解法
（4）数值积分解法
2、所有知识点代码
3、结果---以三阶Runge-Kutta公式为例（其他的类似）
1、概述
（1）常微分方程初值问题数值解法
（2）解题步骤
（3）数值微分解法
（4）数值积分解法
2、所有知识点代码
import numpy as np
import matplotlib.pyplot as pltdef funEval(x,y): #近似值fxy = (x*y-y**2)/x**2 #1#fxy=2*y/x+x**2*np.e**x #2#stablity#fxy=-30*y#fxy= np.e**(-x**2)#fxy=1-y#fxy=2*y/x+x**2*np.e**xreturn fxy
def funtrue(x): #真实值ft=x/(0.5+np.log(x)) #1#ft=x**2*(np.e**x-np.e) #2#stablityu#ft = np.e**(-30*x) #ft = 1-np.e**(-x) #ft=x**2*(np.e**x-np.e)return ft
def Euler(a,b,f,y0,n): #Ｅｕｌｅｒ公式h=np.abs(b-a)/(n-1)y=np.zeros((n,1))x=np.zeros((n,1))y[0] =y0x[0] = afor i in range(1,n,1):x[i]=a+i*hy[i]= y[i-1]+h*f(x[i-1],y[i-1])return x,y
def ModEuler(a,b,f,y0,n): #改进Ｅｕｌｅｒ公式h=np.abs(b-a)/(n-1)y=np.zeros((n,1))x=np.zeros((n,1))y[0] =y0x[0] = afor i in range(1,n,1):x[i]=a+i*hy[i]= y[i-1]+h*f(x[i-1],y[i-1])y[i] = y[i-1]+h/2*(f(x[i-1],y[i-1])+f(x[i],y[i]))return x,ydef Heun(a,b,f,y0,n): #二阶Ｒｕｎｇｅ—Ｋｕｔｔａ方法：Ｈｅｕｎ公式h=np.abs(b-a)/(n-1)y=np.zeros((n,1))x=np.zeros((n,1))y[0] =y0x[0] = aK1,K2=0,0for i in range(1,n,1):x[i]=a+i*hK1 = f(x[i-1],y[i-1])K2 = f(x[i-1]+2/3*h,y[i-1]+2/3*h*K1)y[i] = y[i-1]+h/4*(K1+3*K2)return x,ydef Ord3Kutta(a,b,f,y0,n): #三阶Kutta公式h=np.abs(b-a)/(n-1)y=np.zeros((n,1))x=np.zeros((n,1))y[0] =y0x[0] = aK1,K2,K3=0,0,0for i in range(1,n,1):x[i]=a+i*hK1 = f(x[i-1],y[i-1])K2 = f(x[i-1]+1/2*h,y[i-1]+1/2*h*K1)K3 = f(x[i-1]+h,y[i-1]-h*K1+2*h*K2)y[i] = y[i-1]+h/6*(K1+4*K2+K3)return x,ydef Ord3Heun(a,b,f,y0,n): #三阶Ｈｅｕｎ公式h=np.abs(b-a)/(n-1)y=np.zeros((n,1))x=np.zeros((n,1))y[0] =y0x[0] = aK1,K2,K3=0,0,0for i in range(1,n,1):x[i]=a+i*hK1 = f(x[i-1],y[i-1])K2 = f(x[i-1]+1/3*h,y[i-1]+1/3*h*K1)K3 = f(x[i-1]+2/3*h,y[i-1]+2/3*h*K2)y[i] = y[i-1]+h/4*(K1+3*K2)return x,ydef Ord4Kutta(a,b,f,y0,n): #四阶古典Ｒｕｎｇｅ—Ｋｕｔｔａ公式h=np.abs(b-a)/(n-1)y=np.zeros((n,1))x=np.zeros((n,1))y[0] = y0x[0] = aK1,K2,K3,K4 = 0,0,0,0for i in range(1,n,1):x[i]=a+i*hK1 = f(x[i-1],y[i-1])K2 = f(x[i-1]+1/2*h,y[i-1]+1/2*h*K1)K3 = f(x[i-1]+1/2*h,y[i-1]+1/2*h*K2)K4 = f(x[i-1]+h,y[i-1]+h*K3)y[i] = y[i-1]+h/6*(K1+2*K2+2*K3+K4)return x,ydef Ord4Kutta2(a,b,f,y0,n): #四阶Ｋｕｔｔａ公式h=np.abs(b-a)/(n-1)y=np.zeros((n,1))x=np.zeros((n,1))y[0] = y0x[0] = aK1,K2,K3,K4 = 0,0,0,0for i in range(1,n,1):x[i]=a+i*hK1 = f(x[i-1],y[i-1])K2 = f(x[i-1]+1/3*h,y[i-1]+1/3*h*K1)K3 = f(x[i-1]+2/3*h,y[i-1]-1/3*h*K1+h*K2)K4 = f(x[i-1]+h,y[i-1]+h*K1-h*K2+h*K3)y[i] = y[i-1]+h/8*(K1+3*K2+3*K3+K4)return x,ydef ExpAdamsK1(a,b,f,y0,n): #二步显式Adamsh=np.abs(b-a)/(n-1)y=np.zeros((n,1))x=np.zeros((n,1))y[0] = y0x[0] = afor i in range(1,n,1):x[i]=a+i*hif i ==1:y[i] =y[i-1]+h*f(x[i-1],y[i-1])else:y[i] = y[i-1]+h/2*(3*f(x[i-1],y[i-1])-f(x[i-2],y[i-2])) return x,ydef ExpAdamsK2(a,b,f,y0,n): #三步显式Adamsh=np.abs(b-a)/(n-1)y=np.zeros((n,1))x=np.zeros((n,1))y[0] = y0x[0] = aif n>=2:for i in range(1,n,1):x[i]=a+i*hif i <=2:y[i] =y[i-1]+h*f(x[i-1],y[i-1])else:y[i] = y[i-1]+h/12*(23*f(x[i-1],y[i-1])-16*f(x[i-2],y[i-2])+5*f(x[i-3],y[i-3])) else:print("n must be larger than 2")return x,ydef ModifiedAdamsK2(a,b,f,y0,n): #预测-校正法：四阶Adams预测-校正法h=np.abs(b-a)/(n-1)y=np.zeros((n,1))x=np.zeros((n,1))y[0] = y0x[0] = atemp=0if n>4:x1,y1=Ord4Kutta2(a,h*4,f,y0,4)y[0:4]=y1for i in range(4,n,1):x[i]=a+i*htemp = y[i]+h/24*(55*f(x[i-1],y[i-1])-59*f(x[i-2],y[i-2])+37*f(x[i-31],y[i-3])-9*f(x[i-4],y[i-4]))y[i] = y[i-1]+h/24*(9*temp+19*f(x[i-1],y[i-1])-5*f(x[i-2],y[i-2])+f(x[i-3],y[i-3])) else:print("n must be larger than 4")return x,y
##y'=-30y,y(0)=1, 0<=x<=0.6
def ImplictEuler(a,b,n): #显式Euler公式和隐式Euler法精确度的比较h=np.abs(b-a)/(n-1)y=np.zeros((n,1))x=np.zeros((n,1))y0 = 1y[0] = y0x[0] = aK1,K2,K3,K4 = 0,0,0,0for i in range(1,n,1):x[i]=a+i*hy[i] = y[i-1]/(1+30*h)yt = np.e**(-30*x)return x,y,yt
def main():a,b = 1,1.5 #12n = [2,5,10,20,40,80,160,320,640]#n=[5]y0 = 2 #1#y0 = 0 #2#stablity#a,b = 1,2#y0 = 0# for i in n:# x,y=Euler(a,b,funEval,y0,i)# plt.plot(x,y,label='Euler-solution-'+str(i))# plt.plot(x,funtrue(x),label='True-solution-'+str(i))# plt.legend()# plt.show() for i in n:#x,y=ModEuler(a,b,funEval,y0,i)#x,y=ModifiedAdamsK2(a,b,funEval,y0,i)#x,y,yt=ImplictEuler(0,0.6,i)x,y=Ord3Kutta(a, b, funEval, y0, i)yt= funtrue(x)plt.plot(x,y,label='Euler-solution-'+str(i))plt.plot(x,yt,label='True-solution-'+str(i))plt.legend()plt.show()print(x[0:5],(y-yt)[0:5])print(y)if __name__ == '__main__':main()3、结果---以三阶Runge-Kutta公式为例（其他的类似）
[[1. ][1.5]] [[ 0. ][-0.03207385]]
[[2. ][1.62453333]]
常微分方程初值问题数值解法[完整公式]（Python）相关推荐
- 【数理知识】《数值分析》李庆扬老师-第9章-常微分方程初值问题数值解法
第8章 回到目录 无 第9章-常微分方程初值问题数值解法 9.1 引言 利普希茨 (Lipschitz) 条件 / 利普希茨常数 定理1 解的存在唯一性定理 定理2 解对初值依赖的敏感性 9.2 简单 ...
- 一阶常微分方程的数值解法（二阶显式、隐式 Adams 公式及 Milne 方法）
一阶常微分方程的数值解法 这里我们介绍二阶显式.隐式 Adams 公式及 Milne 方法求解方程. 题目: 对初值问题 u ′ = u − t 2 , 0 ≤ t ≤ 1 , u ( 0 ) = 0 ...
- 常微分方程数值解的c语言程序,常微分方程的数值解法 一阶常微分方程数值解的C语言编程实现.doc...
常微分方程的数值解法 一阶常微分方程数值解的C语言编程实现 导读:就爱阅读网友为您分享以下"一阶常微分方程数值解的C语言编程实现"资讯,希望对您有所帮助,感谢您对92的支持! 一阶 ...
- [常微分方程的数值解法系列三] 改进欧拉法(预估校正法)
改进欧拉法 简介 预估-校正 截断误差 例子 在惯性导航以及VIO等实际问题中利用IMU求解位姿需要对IMU测量值进行积分得到需要的位置和姿态,其中主要就是求解微分方程.但之前求解微分方程的解析方法主 ...
- [常微分方程的数值解法系列二] 欧拉法
欧拉法 简介 几何意义 证明 泰勒展开近似 求导近似 积分近似 几种欧拉方式 向前欧拉公式 向后欧拉公式 梯形公式 中点公式 截断误差 求解过程 向前欧拉公式 例子 向前欧拉公式 在惯性导航以及VIO ...
- [常微分方程的数值解法系列四] 中值法
中值法 简介 具体步骤 截断误差 例子 在惯性导航以及VIO等实际问题中利用IMU求解位姿需要对IMU测量值进行积分得到需要的位置和姿态,其中主要就是求解微分方程.但之前求解微分方程的解析方法主要是应 ...
- 数值分析 第七章 常微分方程的数值解法
1 数值解法相关公式 1.1 为什么要研究数值解法? 所谓数值解法,就是设法将常微分方程离散化,建立差分方程,给出解在一些离散点上的近似值. 1.2 问题 7.1 一阶常微分方程初值问题的一般形式 { ...
- [常微分方程的数值解法系列五] 龙格-库塔(RK4)法
龙格-库塔法 简介 基本思想 具体方法 一阶 二阶 求解参数 特殊二阶 三阶 高阶 步长选择 例子 在惯性导航以及VIO等实际问题中利用IMU求解位姿需要对IMU测量值进行积分得到需要的位置和姿态,其 ...
- 二阶龙格库塔公式推导_[常微分方程的数值解法系列五] 龙格-库塔(RK4)法
在惯性导航以及VIO等实际问题中利用IMU求解位姿需要对IMU测量值进行积分得到需要的位置和姿态,其中主要就是求解微分方程.但之前求解微分方程的解析方法主要是应用于一些简单和特殊的微分方程求解中,对于 ...
最新文章
- WebAssembly能不能取代JavaScript？15张卡通图给你答案！
- 基于VTK的MFC应用程序开发(3)
- 51单片机开发板(W25Q16学习)
- iOS网络请求安全（JWT，RSA）
- Kali渗透(二)之被动信息收集
- 荣耀9igoogle模式_iGoogle个性化主页的6种替代方法
- 大数据Hadoop集群中常用的任务调度框架
- LeetCode 33. Search in Rotated Sorted Array
- Android开发网络连接超时
- MySQL 5.5 使用 Event定期自动维护/执行Procedure
- WebX框架使用说明
- Android 版本检测更新
- 【神兵利器】介绍一款支持屏幕录制、滚动截图、高清长图、图片编辑、图片转PDF格式、屏幕取色的截图软件：FastStone Capture
- crentso7.4+rpm方式安装MySQL5.7.22报错：安装冲突conflicts
- 《拆掉思维里的墙》读后感
- 【地图易-制图案例】全球地震分布地图
- 笔记本服务器管理器在哪个文件夹,笔记本云服务器在哪个文件夹
- 余三码的优点及其与8421码的对比
- marve register license
- 超融合架构和服务器虚拟化是什么关系？主流超融合厂商服务器虚拟化产品对比分析