小刻也能看懂的常微分方程入门!
基本上贺的 Ordinary Differential Equations A Introductiont to the Fundamentals (Kenneth B. Howell) 的内容,看个乐子就好
感谢MMA对其中一些运算的大力支持
基本概念
什么是微分方程,什么是阶数,什么是微分方程的解就略过了
一般要求微分方程解在给定的区间应当是连续的
一些基本应用
物理里有很多微分方程的例子
例如最简单的物体仅收重力的情况下下落位移 y 和时间 t 的关系
根据物理关系 dtdy 是 t 时刻的速度 v(t) dt2d2y 是 t 时刻的加速度 a(t)
显然 a(t)=−g 是恒成立的
假设物体从 y0 处无初速度释放,则 y(0)=y0 同时 v(0)=0
所以 ma(t)=−mg=mdt2d2y
⇔dtd2y=−gdt⇒∫dt2d2ydt=−∫gdt⇒dtdy+C1=−gt+C2
⇒dtdy=−gt+C⇒y(t)=−2gt2+C1x+C2
把两个初值条件带入得到
y(t)=−2gt2+y0
当然实际上空气阻力也是个不可忽略的重要因素,尝试调整方程加入空气阻力项
我们暂且假定空气阻力仅与速度成正比关系 (实际上应该是二次吧)
既 mdtdv=−mg−γv
⇒v(t)=−gt−mγ∫v(t)dt+C
同时利用 v(t) 与 y(t) 的关系可以得到 v(t)=−gt−mγy(t)+C
作为 Basic 部分原书到此为止,不过这个方程还是相当好解的
v(t)=−b(e−bt−1)g(b=mγ)
积分和微分方程
对于形如 dxndny=F(x) 的形式,只要反复积分就行了
当然每次积分都会出现一个任意的系数 Ck 所以需要初值条件来确定具体的的函数
当然用定积分也是可以的,因为 ∫axdxdydx=y(x)−y(a)
一般来说如果给出了 y(x0)=y0 的话,取 a=x0
例如 dxdy=3x2y(2)=12
则 ∫2xdxdydx=∫2x3x2dx
得到 y(x)=x3−8+y(2)=x3+4
用定积分可以避免常数 C 的出现,更重要的是如果无法不定积分出解析解,仍然可以进行计算
例如 dxdy=e−x2y(0)=0
这个显然是积不出来的,但是利用定积分仍然可以计算某一点的值,实际应用上非常重要
因为积不出来很正常,所以定义了一些常用的积分函数
erf(x)=∫0xπ2e−s2ds
Si(x)=∫0xssinsds
因为这种函数性质研究的比较好,所以解出来包含这样的形式一般也是可以接受的
如 dxdy=e−x2y(0)=0 就可以解出来 y(x)=2πerf(x) 也是能接受的
分段函数的微分方程
具体的处理上并不困难,只要分段做就可以,理论上的问题在于对一个分段的函数积分出来的函数是否连续,答案是肯定的
假设函数 f(s) 在 (a,b) 上除了有限的点外均连续,且每个不连续点的左右极限均存在且有限
则积分函数 g(x)=∫axf(s)ds 在 (a,b) 上连续,依照可积性理论不难确定 g(x) 是存在的
所以只需证明 x→x0limg(x)=g(x0) 对于每个间断点 x0 成立即可
根据定义
x→x0limg(x)=x→x0lim∫axf(s)ds=x→x0lim(∫ax0f(s)ds+∫xx0f(s)ds)=g(x0)+x→x0lim∫x0xf(s)ds
所以只需证 x→x0lim∫x0xf(s)ds 即可
根据条件在任意 [α,β]∈(a,b) 均 ∃M,f(x)≤M 成立,所以
0≤x→x0+lim∣∣∣∣∣∫x0xf(s)ds∣∣∣∣∣≤x→x0+limM(x−x0)=0
同理可证左极限
一阶常微分方程
这个国内高数书一般都给出了比较通用的解法,分离变量和常数变易
这一部分探索了一阶常微分方程更深入的一些内容,并且主要运用了积分因子法
常数解
顾名思义就是 y(x)=C 的解,也叫做平衡解,因为如果得到这种解,往往就意味着在各种因素的作用下最终 y 稳定
对于 dxdy=F(x,y) 如果存在有常数解,就意味着存在 C 使得 F(x,C)=0 恒成立
例如函数 dxdy=(y−2x)(y2−9)
常数解就是 y=3 或 y=−3
(个人感觉常数解类似于通解的一个极限,例如考虑 dxdy=2xy2−4xy 有常数解 y=2 或 y=0 和一个通解,而通解的常数趋于无穷时就会出现常数解)
解的存在性与唯一性
一般研究的问题都是从实际上出发,所以都有解且唯一
对于形如 dxdy=F(x,y) 且给出初值 y(x0)=y0 的而微分方程,多元函数 F(x,y) 与 ∂y∂F 若在某一个包含 (x0,y0) 的开区域 R 上连续,则方程的解存在且在某一个包含 x0 的开区间 (a,b) 上唯一
考虑证明它,不妨设 x0=0 ,然后尝试把初值条件变成更好处理的形式
利用之前定积分的思路可以导出 y(x)=∫0xF(s,y(s))ds+y0(x∈(a,b))
带入验证也能发现这个条件确定了初值,所以这两种形式等价
于是我们尝试构造一个函数列来"逼近"解,先口胡一个满足初值的函数,比如 ψ0(x)=y0
然后我们用 ψk+1(x)=y0+∫0xψk(s)ds 来构造下一个函数 (均在 (a,b) 上)
不难发现每个函数都满足初值条件,而如果这个函数列按照我们的期望收敛,记 y(x)=k→+∞limψk(x)
y(x)=y0+k→+∞lim∫0xF(s,ψk(s))ds=y0+∫0xk→+∞limF(s,ψk(s))ds=y0+∫0xF(s,y(s))ds
如果这些取极限是正确的,那么这就是满足初值条件的一个解,这个过程被称作 Picard 方法,函数列 ψ(x) 叫做 Picard 列
严格证明这个内容显然需要一些分析学的内容,以下内容不完全来自原书
Lemma 1
在存在性条件下,存在有正常数 M,B 跨越 0 的闭区间 [a,b] 与有限正实数 Δ
使得在矩形 R1={(x,y)∣a≤x≤b,∣y−y0∣≤Δ} 中 ∣F(x,y)∣≤M,∣∣∣∣∣∂y∂F∣∣∣∣∣≤B
且 0<−aM≤Δ,0<bM≤Δ
同时若有函数 ϕ(x) 在 (a,b) 上满足 ∣ϕ(x)−y0∣≤Δ 则 ψ(x)=y0+∫0xF(s,ϕ(s))ds 在闭区间 [a,b] 上连续且满足 ∣ψ(x)−y0∣≤Δ
Lemma 1 的证明
先随便在 R 中取一个矩形 R0 设 Δ=max{y−y0},y∈R0 矩形 x 轴边界为 a0,b0
因为 R 中 F 与 ∂y∂F 就是连续的,那么在这个闭区域里必然有界
记这个界为 ∣F(x,y)∣≤M,∣∣∣∣∣∂y∂F∣∣∣∣∣≤B
为了满足条件,调整这个矩形,设 ΔX=MΔ
取 a=max(a0,−ΔX),b=min(ΔX,b0)
于是 −aM≤ΔXM=Δ,bM≤ΔXM≤Delta
所以这个 x 轴以 a,b 为界的闭矩形 R1 满足条件
最后,如果 ∣ϕ(x)−y0∣≤Δ 则 (s,ϕ(s))∈R1 即 ∣F(s,ϕ(s))∣≤M
并且显然 F(x,ϕ(x)) 是连续的,所以 ψ(x)=y0+∫0xF(s,ϕ(s))ds 存在
关于 ψ(x) 的连续性只需把表达式代入验证即可
同时 ∣ψ(x)−y0∣=∣∫0xF(s,ϕ(s))ds∣≤M∣x∣≤MΔX=Δ
证毕
Picard 函数列的收敛
因为 ψ0(x)=y0 所以在使用 Lemma 1 时 ∣ψ0(x)−y0∣≤Δ 总是成立的
反复应用 Lemma 1 可以得到 ∣ψk(x)−y0∣≤Δ 成立
为了搞出这个东西的敛散性,我们需要知道 ∣ψk+1(x)−ψk(x)∣ 的界
k=0 时 ∣ψ1(x)−ψ0(x)∣=∣ψ1(x)−y0−(ψ0(x)−y0)∣=∣ψ1(x)−y0∣−∣ψ0(x)−y0∣≤2Δ
考虑 k≥1 时
∣ψk+1(x)−ψk(x)∣=∣∣∣∣∣y0+∫0xF(s,ψk(s))ds−y0−∫0xF(s,ψk−1(s)ds∣∣∣∣∣
∣∣∣∣∣∫0x[F(s,ψk(s))−F(s,ψk−1(s))]ds∣∣∣∣∣≤∫0x∣F(s,ψk(s))−F(s,ψk−1(s))∣ds
且由微分中值定理可得 ∣F(s,ψk(s))−F(s,ψk−1(s))∣≤B∣ψk(s)−ψk−1(s)∣
故而 ∣ψ2(x)−ψ1(x)∣≤B∫0x∣ψ1(x)−ψ0(x)∣ds≤B∫0x2Δds≤2ΔBx
依次递推可以得到
∣ψk+1(x)−ψk(x)∣≤2Δk!(B∣x∣)k,x∈[a,b],k∈N
记无穷和 S(x)=k=0∑∞[ψk+1(x)−ψk(x)]
因为 ∑[ψk+1(x)−ψk(x)]≤∑∣ψk+1(x)−ψk(x)∣
而应用泰勒展开式可得
∑∣ψk+1(x)−ψk(x)∣≤k=0∑∞∣ψk+1(x)−ψk(x)∣=k=0∑∞2Δk!(B∣x∣)k=2ΔeB∣x∣
所以无穷和 S(x) 收敛,所以 k→+∞limψk(x)=ψ0(x)+S(x)
所以 y(x)=k→+∞limψk(x) 是存在的
y(x) 的连续性
根据 Lemma 1 ψk(x) 是连续的,而对于 ∣y(x1)−y(x)∣=k→∞lim∣ψk(x1)−ψk(x)∣≤M∣x1−x∣
所以 y(x) 也必然连续
最后的部分
在之前的讨论中 y(x)=y0+∫0xF(s,y(s))ds 需要建立在函数列 {F(s,ψk(s))} 一致收敛的前提下,所以我们只需要证明该函数列一致收敛即可
∣F(s,y(s))−F(s,ψn(s))∣≤B∣y(s)−ψn(s)∣
∣y(s)−ψn(s)∣=∣∣∣∣∣∣k=0∑∞[ψk+1(x)−ψk(x)]−k=0∑n−1[ψk+1(x)−ψk(x)]∣∣∣∣∣∣=∣∣∣∣∣∣k=n∑∞[ψk+1(x)−ψk(x)]∣∣∣∣∣∣
∣∣∣∣∣∣k=n∑∞[ψk+1(x)−ψk(x)]∣∣∣∣∣∣≤k=n∑∞∣ψk+1(x)−ψk(x)∣≤k=n∑∞2Δk!(B∣x∣)k
换个角标让求和从 0 开始
k=0∑∞2Δ(n+k)!(B∣x∣)n+k≤2Δk=0∑∞n!k!(B∣x∣)n+k≤2Δn!(B∣x∣)nk=0∑∞k!(B∣x∣)k=2Δn!(B∣x∣)neB∣x∣
上式用到了 n!k!≤(n+k)! 很容易验证
同时因为 n→∞limn!(B∣x∣)n=0 所以一定能找到某一个 n 使得 ∣F(s,y(s))−F(s,ψn(s))∣≤B∣y(s)−ψn(s)∣≤δ 在 (a,b) 上成立
所以我们先前交换极限和积分是正确的,至此我们完成了 y(x) 就是解的证明
(如果你不理解一致收敛和交换积分极限的关系,可以尝试证明 ∣ψn(x)−y0−∫F(s,y(s))ds∣→0,并不难证明
以及唯一性
唯一性的证明和存在性差别不大,仅简述思路
惯用的手法,假设除了我们得到的解外还有一个解 y2(x)
根据之前推出的条件,y2(x)=y0+∫0xF(s,y2(s))ds
则 ∣y2−y0∣≤∫0x∣F(s,y2(s))∣ds≤M∣x−0∣≤bM≤Δ
所以 y2(x) 满足 Lemma 1 的初始条件,设 ϕ0(x)=y2(x) ,ψ0(x)=y1(x)
这里 y1(x) 是 Picard 序列得到的解
对函数列 {ϕk(x)} 与 {ψk(x)} 仿照之前的证明,可以得到 ∣ψk+1(x)−ϕk+1(s)∣≤∫0xB∣ψk(s)−ϕk(s)∣ds
之后也能得到 k→+∞lim∣ψk(x)−ϕk(x)∣=0
而因为 y1,y2 均满足 y=y0+∫0xF(s,y(s))ds
所以 ψk(x)=y1,ϕk(x)=y2
于是两个解相同,从而证明了唯一性
区间在哪
先前证明的定理虽然告诉了我们解在某个区间存在且唯一,但是并没有指明区间是什么
其他条件与先前定理相同,而进一步的,如果在带状区域 {(x,y)∣a<x<b,y∈R} 上 F 与 ∂y∂F 均连续且 ∂y∂F 取值仅与 x 有关(实际上只要在这个区域有界即可),则微分方程的解在 (a,b) 上存在且唯一
证明思路仍然类似,不写了
可分离变量的一阶方程
也就是 dxdy=f(x)g(y) 的方程,基础思路就是变成 g(y)1dy=f(x)dx,然后分别积分
这部分高数书里应该都有的(
这个过程结果是对的,但是中间有一些小的问题, dxdy 代表的求导运算,应该不能做出来乘 dx 这个操作
实际上的处理手段应该左边是变成 g(y)1dxdydx,这里的 dx 应理解为微分中的 dx
即 g(y)1y′dx=g(y)1dy,因为 y′dx=dy
所以这个变形过程是对的
常数解
当 g(y)=0 的时候这样做没有任何问题,如果存在 y0 使得 g(y0)=0 ,那根据之前的常数解的内容,y(x)=y0 应当是一个常数解
需要注意的是有的时候常数解和通解可以合并
例如 dxdy=2x(y−5),通解为 y=5±eCex2,常数解为 y=5,而 ±eC 可以取到除 0 外任何值,取 0 时即为常数解,所以最终的解为 y=5+Cex2,C∈R
通用解法
对于 dxdy=f(x)g(y) 的形式
先解出 g(y)=0 所有解作为常数解
随后对 g(y)1dy=f(x)dx 两边积分,然后得到解
存在性与唯一性
考虑含有初值条件的可分离微分方程,不妨设给出初值条件 y(x0)=y0
之前的存在性与唯一性定理要求 F(x,y) 在包含 (x0,y0) 的开区域上连续 (还有他的偏微分)
而如果不存在这个开区间意味着可能不止一个解
例如考虑 dxdy=2y,y(1)=0
F(x,y)=2y,∂y∂F=y1 在任意包含 y0=0 的开区间均不连续,所以这个方程有不止一个解
例如常数解 y=0 分段解 y=0 在 x<1 上y=(x−1)2 在 x≥1 上
而且在解满足唯一性的具有初值条件的微分方程具体表达式上也可能会出现其他的解,但是这些解是错误的,例如 dxdy=2y,y(0)=4
它显然仅有一个解,但是在解的过程中,我们会遇到 2y=2x+C 的情况
如果你在这一步带入初值条件解出 C=4 是正确的,但如果你先平方后带入会得到 C=±4 的解
原因是经典的的 A2=∣A∣
总之满足唯一性条件的方程一定唯一,如果解出来多个解一定有某个步骤错了
线性一阶方程
线性一阶方程指 dxdy+p(x)y=f(x) 的微分方程
例如 x2dxdy+x3[y−sin(x)]=0
可以化为 dxdy+xy=xsin(x),即线性方程
解法
考虑不退化的线性方程 dxdy+py=f
假设我们现在有一个函数 μ(x) 满足 dxdμ=μp,则 dxd[μy]=μdxdy+ydxdμ=μdxdy+μpy
在原方程两边同时乘 μ 则 μdxdy+μpy=μf 即 dxd[μy]=μf
我们就可以两边积分得到解了,而 μ 只需要解一个可分离变量的微分方程即可得到 μ=±e∫p(x)dx ,这个函数被称为积分因子,而我们只需要一个满足条件的 μ 即可,所以正负号任取其一
小题狂练(?
xdxdy+4y=x3
不难得到 p(x)=x4,f(x)=x2 则 μ=e4ln∣x∣=∣x4∣=x4
随后 μy=∫μfdx 得到 x4y=71x7+C
最终 y=71x3+Cx−4
另一个小技巧是去掉绝对值,很多时候 μ 计算出来的结果是某个连续的简单函数的绝对值,为了避免绝对值的麻烦,可以直接去掉绝对值(什么是简单函数依赖你的判断力)
试试看,解 dxdy+cot(x)y=xcsc(x)
答案是 y=2sin(x)x2+C
使用积分因子的方法的时候我们仍然可以使用定积分简化初值问题的运算
例如 dxdy−2xy=4,y(0)=3
积分因子 μ=e−x2+C 同样的,这个 C 无关紧要,取 0 来简化运算
之后
y(x)=μ1[μ(0)y(0)+∫0xμ(s)f(s)ds]
最终为 y(x)=ex2(3+4∫0xe−s2ds)=ex2[3+2πerf(x)]
存在性与唯一性
重新回到这个话题,这里给出的方法实际上和原方程等价,所以给定初值情况下有且仅有这一个解
更宽松一点,只要函数 p,f 有这有限多个间断点,方程解的存在和唯一都不受影响
间章:换元法
在进入到更一般的微分方程之前,我们先看看我们能力的限度
对于可以直接积分,可以分离变量和线性一阶方程我们都可以简单的计算出一个满意的结果
但是仍然有相当多的问题我们没有成体系的方法
例如 dxdy=(x+y)2 这个简单的方程不适用任何我们掌握的方法
而将 y 替换成合适的函数 可能会带来些许帮助
设 u=x+y 则 y=u−x 推出 dxdy=dxdu−dxdx=dxdu−1
即 dxdu=u2+1 这个就很容易解出来 u=tan(x+C)
所以 y=tan(x+C)−x
有些换元没啥套路,很难搞,但是有一些还是能搞定的
线性换元
对于 dxdy=f(Ax+By+C) 的形式,尝试 u=Ax+By+C 则可以得到 dxdu=A+Bf(u)
例题 dxdy=2x−4y+71
答案是 y=2x+45 (来自 u 的常数解)或 ln∣2x−4y+5∣=4y+C
比值换元
对于 dxdy=f(xy) 的形式,利用 u=xy 代换,得到 dxdu=xf(u)−u
例题 xy2dxdy=x3+y3
答案为 y=x33ln∣x∣+C
伯努利等式
对于 dxdy+p(x)y=f(x)yn
一般的直觉应该是用 u=yn 换元,但是稍加思索会发现这样行不通,这样不会简化任何东西
而我们如果用 u=yr 换元,化简整理能得到 ru1+rn−1f=rpu+dxdu
只要取 r=1−n 就能得到 (1−n)f=(1−n)pu+dxdu 的一阶线性微分形式
而显然 y=0 也是一个解
例题 dxdy+6y=30e3xy2/3
答案为 y=0 与 y=[2e3x+Ce−2x]3
(答案为啥没过程?因为用MMA算的)
微分方程的恰当形式与一般解法
微分方程的恰当形式( Exact Form )是指对于微分方程 M(x,y)+N(x,y)dxdy=0
存在可微函数 ϕ(x,y) 满足 ∂x∂ϕ=M(x,y),∂y∂ϕ=N(x,y) 恒成立
借由向量场中的理论 ϕ 被叫做势函数
显然不是所有微分方程都由恰当形式的,而由恰当形式的微分方程也不一定仅有一个势函数
但是由之前的理论,给定初值后解的存在应当在一定条件下唯一,所以选取什么势函数并不重要
解具有恰当形式的微分方程
假设我们已知某个微分方程的一个势函数 ϕ(x,y)
则微分方程可以写为 ∂x∂ϕ+∂y∂ϕdxdy=0 逆用链式法则,得到
∂x∂ϕ(x,y(x))=0,两边积分得到 ϕ(x,y(x))=C 也就是 ϕ(x,y(x))=C 就是原方程的解
例如 2xy+(2y+x2)dxdy=0
一个势函数 ϕ(x,y)=y2+x2y ,因此 y2+x2y=C 就是原方程的一个解
寻找势函数
可是一般不会有小鸟飞进来告诉我们势函数是什么,我们需要建立一个找到势函数的方法
能发现势函数由两个偏微分方程确定,即 ∂x∂ϕ=M(x,y),∂y∂ϕ=N(x,y)
我们用例子 2xy+2+(x2+4)dxdy=0 来说明这个过程
先看 ∂x∂ϕ=2xy+2 这个方程两边对 x 积分,得到 ϕ(x,y)=x2y+2x+h(y)
现在我们只要知道 h(y) 即可,将现在的 ϕ(x,y) 带入 ∂y∂ϕ=N(x,y) 中
得到 x2+h′(y)=x2+4 导出 h′(y)=4 现在我们能轻松解出来 h(y)=4y+C
需要注意的是,如果在这里没办法得到一个只含 y 的方程,那就意味着这个方程并不是恰当形式
势函数的存在性
寻找势函数的过程中虽然能够判断势函数是否存在,但是这样的判断往往需要耗费不少力气
而由微分学 ∂x∂y∂2ϕ=∂y∂x∂2ϕ 所以 ∂y∂M=∂x∂N 应当成立,这是一个必要条件
当这个条件成立时就可以试试寻找势函数(实际上这是一个充要条件,这个条件看起来应该有点眼熟)
将非恰当形式转变为恰当形式
考虑方程 3y+3y3+(xy2−x)dxdy=0
可以判断这个方程并不具有恰当形式,当时如果我们乘上一个 x2y−2 就能变成
3x2y−1+3x2y+(x3−x3y−2)dxdy=0
而这就是一个恰当形式
考虑我们乘的这个东西(不妨叫他恰当因子?原书称他为一般积分因子)
我们的目标是找到 μ(x,y) 使得 ∂y∂μM=∂x∂μN
这个讨论超出了目前话题范围,所以仅考虑一些特殊的积分因子
∂y∂μM=∂x∂μN⇒μ∂y∂M+M∂y∂μ=μ∂x∂N+N∂x∂μ
⇒μ(∂y∂M−∂x∂N)=N∂x∂μ−M∂y∂μ
考虑如果 ∂y∂μ=0 意味着 μ 仅与 x 有关,此时 μ1dxdμ=N∂y∂M−∂x∂N
记右边为 ϕ(x) 则 μ=e∫ϕ(x)dx 即一个积分因子,而仅与 y 有关是结论类似此时 ϕ(y)=−M∂y∂M−∂x∂N
如果不是这两种情况最好试试 xayb 的形式
而如果 μ 是 函数 δ(x,y) 的一个函数,则
lnμ(δ)=∫M∂y∂δ−N∂x∂δ∂x∂N−∂y∂Mdδ
当然找的这个 δ 则是另一个困难的问题
至此对于一阶微分方程性质上的讨论基本完结,而实际应用上更多的需要一种数值算法来迭代逼近解,这部分内容留在之后的 1.5 讨论(如果有的话