常微分方程初值问题数值解法[完整公式](Python)
Python 完整 数值 公式 解法 微分方程
2023-09-14 09:05:25 时间
目录
3、结果---以三阶Runge-Kutta公式为例(其他的类似)
1、概述
(1)常微分方程初值问题数值解法
(2)解题步骤
(3)数值微分解法
(4)数值积分解法
2、所有知识点代码
import numpy as np
import matplotlib.pyplot as plt
def 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**x
return 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): #Euler公式
h=np.abs(b-a)/(n-1)
y=np.zeros((n,1))
x=np.zeros((n,1))
y[0] =y0
x[0] = a
for i in range(1,n,1):
x[i]=a+i*h
y[i]= y[i-1]+h*f(x[i-1],y[i-1])
return x,y
def ModEuler(a,b,f,y0,n): #改进Euler公式
h=np.abs(b-a)/(n-1)
y=np.zeros((n,1))
x=np.zeros((n,1))
y[0] =y0
x[0] = a
for i in range(1,n,1):
x[i]=a+i*h
y[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,y
def Heun(a,b,f,y0,n): #二阶Runge—Kutta方法:Heun公式
h=np.abs(b-a)/(n-1)
y=np.zeros((n,1))
x=np.zeros((n,1))
y[0] =y0
x[0] = a
K1,K2=0,0
for i in range(1,n,1):
x[i]=a+i*h
K1 = 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,y
def 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] =y0
x[0] = a
K1,K2,K3=0,0,0
for i in range(1,n,1):
x[i]=a+i*h
K1 = 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,y
def Ord3Heun(a,b,f,y0,n): #三阶Heun公式
h=np.abs(b-a)/(n-1)
y=np.zeros((n,1))
x=np.zeros((n,1))
y[0] =y0
x[0] = a
K1,K2,K3=0,0,0
for i in range(1,n,1):
x[i]=a+i*h
K1 = 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,y
def Ord4Kutta(a,b,f,y0,n): #四阶古典Runge—Kutta公式
h=np.abs(b-a)/(n-1)
y=np.zeros((n,1))
x=np.zeros((n,1))
y[0] = y0
x[0] = a
K1,K2,K3,K4 = 0,0,0,0
for i in range(1,n,1):
x[i]=a+i*h
K1 = 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,y
def Ord4Kutta2(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] = y0
x[0] = a
K1,K2,K3,K4 = 0,0,0,0
for i in range(1,n,1):
x[i]=a+i*h
K1 = 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,y
def ExpAdamsK1(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] = y0
x[0] = a
for i in range(1,n,1):
x[i]=a+i*h
if 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,y
def ExpAdamsK2(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] = y0
x[0] = a
if n>=2:
for i in range(1,n,1):
x[i]=a+i*h
if 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,y
def 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] = y0
x[0] = a
temp=0
if n>4:
x1,y1=Ord4Kutta2(a,h*4,f,y0,4)
y[0:4]=y1
for i in range(4,n,1):
x[i]=a+i*h
temp = 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 = 1
y[0] = y0
x[0] = a
K1,K2,K3,K4 = 0,0,0,0
for i in range(1,n,1):
x[i]=a+i*h
y[i] = y[i-1]/(1+30*h)
yt = np.e**(-30*x)
return x,y,yt
def main():
a,b = 1,1.5 #12
n = [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 re.compile() 详解——Python正则表达式「建议收藏」
- pycharm创建python虚拟环境好处_pycharm虚拟环境与本地环境区别
- 用Python做一个游戏辅助脚本,完整编程思路分享
- python 生成数组_Python创建数组「建议收藏」
- python移动app开发_神奇的Kivy,让Python快速开发移动app
- Python入门系列(五)一篇搞懂python语句
- Python入门系列(六)一篇学会python函数
- doc转docx和docx转pdf的python代码2021.9.28
- OpenCV-Python学习(15)—— OpenCV 鼠标操作和响应(cv.setMouseCallback)
- 【说站】python三位数反序输出
- 【说站】python读取txt文件
- 【说站】python tempfile创建文件
- python安装numpy后pycharm导入不了_如何导入numpy
- python 类属性和实例属性、类方法, 静态方法, 实例方法、接口, 协议和抽象基类 (4.2)
- python pkl文件_Python字符串格式化输出的方式包括
- h5 Python_python做h5网站
- python lambda表达式举例_Python中lambda表达式[通俗易懂]
- 【python量化】用python搭建一个股票舆情分析系统
- Python绘制旭日图_python绘制散点图
- python如何生成随机数_Python生成50个随机数
- 10 个杀手级的 Python 自动化脚本
- python 按行数分割文件
- Python 自动化指南(繁琐工作自动化)第二版:一、PYTHON 基础知识
- Python unittest模块实现单元测试
- 从Python连接Oracle数据库介绍(python连接oracle)
- 在Linux上运行Python脚本的简单指南(linux运行python)
- Python如何使用MySQL构建立连接(python怎么连接mysql)
- python提取文件的小程序
- Python中apply函数的用法实例教程
- python字符串排序方法