高等数学--微分方程的求解(sympy库)

微分方程的意义

函数是客观事物内部联系的量化,利用函数关系可以研究客观事物的内在规律,因此,如何寻找函数关系在实践中有着极其重要的意义!
但是在很多问题当中,往往是若干“因”导致若干“果”,想要通过观察直接找出变量之间的函数关系比较困难,那么我们该怎么办呢?
我们可以通过问题中反映的情况,列出含有要找的函数及其导数关系式,这样的关系式就是微分方程
微分方程建立之后,需要对其进行研究进而找到未知函数来,这一过程就是解微分方程

引自《高等数学 第七版 同济大学出版社》

sympy求解微分方程的方式

symbols()+Eq()

求解下面微分方程
y ′ ′ − 2 y ′ ′ + y = 0 y''-2y''+y=0 y′′2y′′+y=0

import sympy as sp

x=sp.symbols('x')
y=sp.symbols('y',cls=sp.Function)

eq1 = sp.Eq(y(x).diff(x, x) - 2*y(x).diff(x) + y(x), 0)
res1=sp.dsolve(eq1,y(x))
print(res1)
>>> Eq(y(x), (C1 + C2*x)*exp(x))

Function()

求解下面微分方程
y ′ ′ − 2 y ′ + y = 0 , 其中 y ( 0 ) = 0 , y ′ ( 0 ) = 1 y''-2y'+y=0,其中y(0)=0,y'(0)=1 y′′2y+y=0,其中y(0)=0y(0)=1

from sympy.abc import x
import sympy as sp
y=sp.Function('y')
eq=sp.diff(y(x),x,2)-2*sp.diff(y(x),x)+y(x)

con={y(0):0,sp.diff(y(x),x).subs(x,0):1}
res=sp.dsolve(eq,ics=con)
res
>>> Eq(y(x), x*exp(x))

solveset()

求解下面微分方程
y ′ + y + y 2 = 0 ,其中 y ( 0 ) = 1 y'+y+y^2=0,其中y(0)=1 y+y+y2=0,其中y(0)=1

from sympy import *
f = symbols('f', cls=Function)
x = symbols('x')
eq = Eq(f(x).diff(x,1)+f(x)+f(x)**2, 0)
dsolve(eq, f(x))
>>>  Eq(f(x), -C1/(C1 – exp(x)))
C1 = symbols('C1')
eq1 = -C1/(C1 - exp(x))
eq2 = eq1.subs(x, 0)
solveset(eq2 - 1, C1)
>>> {1/2}

求解实际问题

知道了如何使用sympy库求解各种各样的微分方程,下面我们来解决一点现实中遇到的问题

1. 降落伞问题

某降落伞在下降过程中,受到的空气阻力与速度成正比,设降落伞的初始速度为0,求降落伞下降速度与时间的函数关系

假设变量

降落伞的速度:v(t)
重力:P,大小为mg,方向与v一致,
阻力:R,大小为kv,方向与v相反
降落伞所受到的合力为:F

建立模型

根据题目中给出的比那变量之间的函数关系可得
F = m g − k v F=mg-kv F=mgkv
根据牛顿第二定律
F = m a , t 时刻降落伞的加速度 a = d s 2 d 2 t = d v d t F=ma,\\t时刻降落伞的加速度a=\frac{ds^2}{d^2t}=\frac{dv}{dt} F=mat时刻降落伞的加速度a=d2tds2=dtdv

那么列出微分方程:
m d v d t = m g − k v ,其中 v t = 0 = 0 m\frac{dv}{dt}=mg-kv,其中v_{t=0}=0 mdtdv=mgkv,其中vt=0=0

求解微分方程模型

d v m g − k v = d t m \frac{dv}{mg-kv}=\frac{dt}{m} mgkvdv=mdt
两端积分得到
∫ 1 m g − k v d v = ∫ 1 m d t \int \frac{1}{mg-kv}dv=\int \frac{1}{m}dt mgkv1dv=m1dt

由于mg-kv>0,则
− 1 k l n ( m g − k v ) = t m + C 1 -\frac{1}{k}ln(mg-kv)=\frac{t}{m}+C_1 k1ln(mgkv)=mt+C1


m g − k v = e − k m t − k C 1 mg-kv=e^{-\frac{k}{m}t-kC_1} mgkv=emktkC1

于是
v = m g k ( 1 − e − k m t ) v=\frac{mg}{k}(1-e^{-\frac{k}{m}t}) v=kmg(1emkt)

sympy求解

from sympy import *

m,g,k,t=symbols('m g k t')
v=Function('v')
eq=Eq(m*diff(v(t),t),m*g-k*v(t))
con={v(0):0}
res=dsolve(eq,ics=con)
res
>>> Eq(v(t), g*m/k - g*m*exp(-k*t/m)/k)

从上面我们可以看出,我们在求解微分方程时经过这么多步骤推导出的 v 与 t v与t vt的函数关系式,被sympy的几行代码就解决了!不得不叹服sympy库的便利

2. 放射性元素衰变问题

镭元素的衰变规律如下,衰变速度与其现存量R成正比现在已知1600年后,镭的衰变变为原来质量的 1 2 \frac{1}{2} 21,求镭元素现存量R与时间t之间的函数关系

假设变量

现存量R
衰变常数 λ \lambda λ
衰变速度v

构建与求解模型

根据题目中的已知关系

v = d R d t = − λ R v=\frac{dR}{dt}=-\lambda R v=dtdR=λR
那么
d R R = − λ d t \frac{dR}{R}=-\lambda dt RdR=λdt
两边积分得到
∫ 1 R d R = ∫ − λ d t \int \frac{1}{R}dR=\int -\lambda dt R1dR=λdt

l n R = − λ t + l n C   R = C e − λ t lnR=-\lambda t+lnC\\ \ R=Ce^{-\lambda t} lnR=λt+lnC R=Ceλt
由于
t = 0 时, R = R 0 t = 1600 , R = 1 2 R 0 t=0时,R=R_0 \\t=1600,R=\frac{1}{2}R_0 t=0时,R=R0t=1600,R=21R0
代入上式可得
λ = l n 2 1600 = 0.000433 \lambda=\frac{ln2}{1600}=0.000433 λ=1600ln2=0.000433

于是
R = R 0 e − 0.000433 t R=R_0e^{-0.000433t} R=R0e0.000433t

sympy求解

from sympy import *
t, k ,R0= symbols('t,k,R0')
R = Function('R')
R_ = Derivative(R(t),t)
Eqn = Eq(R_,-k*R(t))
res = dsolve(Eqn,ics={R(0):R0})
print(res)
>>> Eq(R(t), R0*exp(-k*t))
kk=solve(res,k)
k0=kk[0].subs({t:1600,R(t):0.5*R0})# 求式子中的参数k(元素衰变系数)
print(kk,k0)
>>> [log(R0/R(t))/t] 0.000433216987849966
k1=res.args[1]
k11=k1.subs({t:99999,k:k0})
print(k11)
>>> 1.53395780934154e-19*R0   # 99999年之后元素的剩余量
  • 5
    点赞
  • 25
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 2
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

夺笋123

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值