python解常微分方程(组)
本人目前初三,能力所限,如有不足之处,还望多多指教。
一周前看到了一个视频,于是我便想用python来求解这个问题。
〇、分析
假设在平面内有一带电粒子q,质量为m。空间内存在匀强磁场B,方向垂直于平面向内即沿z轴负半轴,以及一个沿y轴负半轴的重力场。带电粒子从磁场内O点释放。
则可直接列出粒子的运动方程
将这个方程分解成x和y两个方向
联立即可求得该方程组的解。
一、sympy中的dsolve方法
#导入
from sympy import *
import numpy as np
import matplotlib.pyplot as plt
init_printing()
首先声明符号x,y,q,m,B,g
#参量
q,m,B,t,g = symbols('q m B t g',real=True,nonzero=True,positive=True)
#函数
x,y = symbols('x y',real=True,nonzero=True,positive=True,cls=Function)
再将微分方程表示出来
eq1 = Eq(x(t).diff(t,t),-q*B*y(t).diff(t)/m)
eq2 = Eq(y(t).diff(t,t),-g+q*B*x(t).diff(t)/m)
sol = dsolve([eq1,eq2])
现在打印出sol(用jupyter notebook效果更好)
sol
很明显,这个式子非常冗杂,用trigsimp()方法化简
x = trigsimp(sol[0].rhs)
y = trigsimp(sol[1].rhs)
x,y
然后,可以将里面的积分常数算出
#定义积分变量,避免报错
var('C1 C2 C3 C4')
#x(0)=0
e1 = Eq(x.subs({t:0}),0)
#x'(0)=0,要将subs放在diff后面
e2 = Eq(x.diff(t).subs({t:0}),0)
#y(0)=0
e3 = Eq(y.subs({t:0}),0)
#y'(0)=0
e4 = Eq(y.diff(t).subs({t:0}),0)
l = solve([e1,e2,e3,e4],[C1,C2,C3,C4])
l
紧接着,将积分常量替换到x和y里面,我们就得到了解的最终形式
x = x.subs(l)
y = y.subs(l)
x,y
当然,这个解还可以写成这种形式
用plt画图来看看
ts = np.linspace(0,15,1000)
consts = {q:1,B:1,g:9.8,m:1}
fx = lambdify(t,x.subs(consts),'numpy')
fy = lambdify(t,y.subs(consts),'numpy')
plt.plot(fx(ts),fy(ts))
plt.grid()
plt.show()
但是sympy有一个缺点,当微分方程很复杂时,它会直接罢工。
于是另一个新的方法替我们解决了这个问题
二、scipy.integrate中的odeint方法
#导入
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from math import e
我们首先要创建一个函数,使它可以表示这个微分方程。
对于这个函数应该具备怎样的形式,先从一阶微分方程开始
比如说,我们要求解下面这个方程
它的解通过一些简单的分离变量即得
下面用odeint来解出它,令y(1)=1
#一阶
def f(x,y):#将一阶导数用y和x构成的函数来表示dydx = -y/xreturn dydx
#初始条件
y0 = 1
#初值条件是在自变量x的下限取到。如y(1)就要生成一个从1开始的范围
x = np.linspace(1,5,100)
plt.xlim(0,5)
plt.ylim(0,5)
#odeint()方法,tfirst属性指的是函数f的第一个参数是自变量
sol = odeint(f,y0,x,tfirst=True)
plt.plot(x,sol[:,0])
plt.grid()
plt.show()
同理,对于下面的一阶微分方程组
此时我们设向量u=(x,y) (列向量),方程组化为
其实对于任意的一阶微分方程,都可以写成这样的形式
f的含义就是
此时微分方程就可以化为
而f(t,u)正是我们要找的那个函数
odeint中支持向量输入,因此,可以这样构造这个函数
def f(t,u):#x1,x2...xn = ux,y = u#dx1dt=...,dx2dt=...,dxndt=...dxdt = 3*x-x*ydydt = 2*x-y+e**t#dudt = [dx1dt,dx2dt,...dxndt]dudt = [dxdt,dydt]return dudt
#x1(0)=...,x2(0)=...,xn(0)=...
u0=[0,0]
t = np.linspace(0,10,100)
sol = odeint(f,u0,t,tfirst=True)
因此,对于二阶常微分方程
我们首先把这个方程化解成这样的形式
此时这是一个关于dy/dx和y的一个微分方程组
令u=(y,dy/dx),我们有
def f(x,u):y,dydx = ud2ydx2 = np.exp(x)-4*dydx-4*ydudx = [dydx,d2ydx2]return dudx
u0=[0,0]
x = np.linspace(0,10,100)
sol = odeint(f,u0,x,tfirst=True)
对于更高阶的微分方程也是同样处理
此时我们在回过头来看最初的问题,可以轻松的得到
def f(t,r,k,g):
#k = B*q/mx,y,dxdt,dydt = rd2xdt2 = -k*dydtd2tdt2 = -g + k*dxdtreturn [dxdt,dydt,d2xdt2,d2ydt2]
定义一些常量
t = np.linspace(0,15,1000)
k = 1
g = 9.8
使用odeint方法
r0 = [0,0,0,0]
sol = odeint(f,r0,t,tfirst = True,args=(k,g))
将图像画出
plt.plot(sol[:,0],sol[:,1])
plt.grid()
plt.show()
odeint解这个方程也十分精确,与sympy相差无几。
最后我们还可以让他实现动画效果
def update_points(num):point_ani.set_data(s[:,0][num], s[:,1][num])return point_ani,
plt.xlabel('X')
plt.ylabel('Y')
fig,ax = plt.subplots()
plt.plot(s[:,0],s[:,1])
point_ani, = plt.plot(s[:,0][0], s[:,1][0], "ro")
plt.grid(ls="--")
# 生成动画
ani = animation.FuncAnimation(fig, update_points,range(1, len(t)), interval=5, blit=True)
plt.show()
如果要深度学习这两个库的话,还是推荐官方文档
sympy:https://docs.sympy.org/latest/tutorial/index.html
scipy:https://docs.scipy.org/doc/scipy/reference/tutorial/index.html
最近还准备在更新几篇(如果期末不挂的话)
注:我在b站上也发布了这篇文章
python解常微分方程(组)相关推荐
- 用python实现解常微分方程组的简单示例以及用odeint解常微分方程的范例
背景: 包括两个部分,一个是演示怎么自己写代码解常微分方程,另一部分就是示范python怎么调用解常微分方程的函数. 下面的方程组给出洛仑兹引子在三个方向上的速度,求解运动轨迹 软件: python3 ...
- python解常微分方程龙格库_求解常微分方程组初值问题的龙格库塔法分析及其C代码...
求解常微分方程组初值问题的 龙格库塔法分析及其 C 代码 1 .概 述 由高等数学的知识可知,一些特殊类型的常微分方程(组)能够求出给定初 始值的解析解, 而在科学与工程问题中遇到的常微分方程 (组) ...
- python解常微分方程
一.sympy.dsolve 首先,感觉最科学的是用sympy的dsolve解常微分方程,直接贴代码 import sympy as sydef differential_equation(x,f): ...
- python解常微分方程_Python-sympy.dsolve求解常微分方程(组)
这里分别介绍怎么利用sympy.dsolve求解常微分方程和常微分方程组. #首先利用sympy.dsolve求解单个的常微分方程: #代码 from sympy import Function, d ...
- python解常微分方程龙格库_数值常微分方程-欧拉法与龙格-库塔法
大三时候在跳蚤市场闲逛,从一位数学院的学长那里买了一些闲书,最近翻出来刚好有李荣华.刘播老师的<微分方程数值解法>和王仁宏老师的<数值逼近>,结合周善贵老师的<计算物理& ...
- matlab中常微分方法,MATLAB解常微分方程组的解法(好东西要共享)
1:问题 常微分方程的初值问题的标准数学表述为:y'=f(t,y),a<=t<=b,y(a)=y(0) :我们要求解的任何高阶常微分方程都可以用替换法化为上式所示的一阶形式,其中y为向量, ...
- python解常微分方程龙格库_求解二阶常微分方程的RungeKutta四阶方法
我试着做一个简谐振子的例子,它将用龙格-库塔四阶法求解.要求解的二阶常微分方程(ODE)和初始条件为: y''+y=0 y(0)=0和y'(0)=1/pi 范围在0到1之间,共有100步.我用u作为辅 ...
- python解常微分方程龙格库_excel实现四阶龙格库塔法runge-kutta解二阶常微分方程范例.xls...
excel实现四阶龙格库塔法runge-kutta解二阶常微分方程范例,rungekutta,四阶rungekutta法,rungekuttamatlab,四阶rungekutta,rungekutt ...
- python解不等式组_LintCode Python解法
3.统计数字(Digit Count) 计算数字 k 在 0 到 n 中的出现的次数,k 可能是 0~9 的一个值. 首先是,惯用思维,2个循环解决,这样做的时间复杂度为O(n*2) 1 classS ...
最新文章
- 百度网盘的速度又又又又又又被黑了...侮辱性极强...
- 报文在三次握手过程中丢失怎么办?
- Java中ArrayList和LinkedList区别 时间复杂度 与空间复杂度
- innerText,outerText,innerHTML,outerHTML区别
- 河北2018年职称计算机开始,2018河北职称计算机考试操作题答案(8页)-原创力文档...
- php提示是否运行,php运行错误提示
- C/C++轻松破解别人程序的窗口标题
- commons-io实现流的拷贝
- 从入门到精通:卷积神经网络初学者指南
- iOS开发问题之Could not instantiate class named NSLayoutConstraint
- vue 实现计算器功能
- java替换图片文字_Java 替换PPT文档中的文本和图片
- 职业学校计算机和机电哪个好,职业学校都有什么专业10大热门专业
- android 临时文件存储,缓存和临时文件/文件夹删除android
- TI MSP430工程配置及2019年电赛A题编程示例(使用430 F5529)
- 如何做好PPT——画图篇
- Shiro限制登录尝试次数(适用于单节点)
- wow语音服务器卡蓝条,魔兽世界6.0卡蓝条解决方法 登陆界面卡主解决办法
- 微信小程序怎么集成腾讯IM
- Sers微服务2.1.1