使用 Julia 求解一阶微分方程

Solving First Order Differential Equations with Julia

Source | HN Comments

文章介绍了使用 Julia 语言的 `DifferentialEquations.jl` 包求解一阶微分方程的方法。文章首先解释了微分方程的概念，然后通过放射性衰变和牛顿冷却定律两个例子，详细演示了如何将一阶微分方程转化为 Julia 代码，并使用该库求解。文章侧重于高级 API 的使用，强调了将数学符号转换为 Julia 可理解形式的过程，并提供了代码示例和结果图示。

Posts by Rito Ghosh Search

使用 Julia 求解一阶微分方程

使用 DifferentialEquations.jl 求解线性一阶 DEs math julia computational-math Author Ritobrata Ghosh Published March 3, 2025

使用 Julia 求解一阶微分方程

Image by Pete Linforth from Pixabay

介绍

在本教程中，我们将学习如何使用 Julia 求解微分方程。微分方程在科学和工程问题中无处不在。使用计算机求解微分方程非常方便。我们将使用 DifferentialEquations.jl 包来求解微分方程。

适用人群

本教程面向已经熟悉数学中的微分方程并且可以编写代码（最好是 Julia 代码）的人。只需要 Julia 语言的基础知识。如果您不熟悉 Julia，但拥有使用 Python、C/C++、Ruby、JavaScript、golang 等任何一种语言的经验，那么您可以轻松掌握本教程。

我们将要做什么

在本教程中，我们将只关注使用 DifferentialEquations.jl 包的高级 API，并且我们将从零开始学习如何将一阶微分方程转换为 Julia 代码，以及如何使用该库来求解它们。我们将不会学习和编写求解微分方程的数值方法，例如 Euler’s method, Runge-Kutta methods 等。我们也不会讨论测试精度、测试算法收敛性的问题。我们仅仅是将一些熟悉的数学符号表示的 DE 转换为 Julia 和库可以理解的形式，然后调用 DifferentialEquations.jl 库的高级 API 来求解我们的 DE。就是这样。在本篇文章中，我们将自己限制在一阶微分方程的范围内。

什么是微分方程：一个非常快速的回顾

“普通”方程

让我们先退一步，讨论一下普通方程。普通方程是包含一个或多个未知根的方程。像这样： x2−3x−18=0 由于这是一个二次方程，因此您可以求解 x 的两个值。你如何_创建_这个方程？您可以从真实世界的场景中构建一个方程。这是一个例子，您将从中得到上面的方程：

一位农民正在设计一个矩形菜园。花园的长度比宽度长 3 米。花园的总面积为 18 平方米。这个花园的边长是多少？这就是方程的创建方式：x(x+3)=18，其中 x 是花园较小边的长度。

微分方程

在普通方程中，我们求解方程以找出我们以前未知的变量值。在微分方程的情况下，我们求解以找出我们以前未知的函数。对于普通方程，我们找出未知变量与我们已知的事物之间的关系，具体地说。在 DE 中，我们知道一个未知函数如何随时间或其他变量_变化_。我们用一个未知函数的导数及其导数（一个或多个）来形成我们的方程。它可能看起来像这样： dydx=−y 微分方程用于对现实世界的事件进行建模，我们的目标是求解该函数。只有当一个函数随另一个变量（例如时间）变化时，才有可能形成微分方程。函数的输出随时间变化是不够的，函数本身也必须随时间变化。这里有一些例子：

一根理想的蜡烛以恒定的速率燃烧，即当蜡烛处于其全部长度时，以及当它几乎完成时，蜡的燃烧速率是相同的。但是，当人口较少时，人口的增长率与人口较高时不同。不是人口本身随着人口的增加而增加，而是人口增长率也随着人口的增长而增加。（这里我说的是现代的人口——在没有天然捕食者，医疗设施非常好的情况下。）
我能想到的另一个例子是用于排空水箱的孔的排水速率。当水箱几乎满了时，从孔中排水的速率非常高，但是当水箱几乎空了时，排水的速率要低得多。
当你刚泡好咖啡，假设它的温度是 70∘C，而你周围的温度是 8∘C 时，咖啡从 70∘C 达到 30∘C 所需的时间比它从 12∘C 达到 8∘C 所需的时间要短得多。当温差较小时，与温差较大时相比，_冷却过程_较慢。

这些是您使用微分方程建模的情况。

微分方程的数值解

在某些情况下，找到微分方程的精确解是不可能的，或者是不切实际的。这就是为什么我们尝试寻找近似解。在许多情况下，这已经足够好了。当无法获得解析解时，数值解是我们所拥有的全部。

我们在这里期望什么

对于我们在这里将要实现并将进行数值求解的微分方程，我们将看到一个图，该图将预测未知函数的行为。

第一个例子：放射性衰变

放射性是其中不稳定的原子核通过辐射失去能量的现象。在放射性衰变中，如果保存一份放射性原子核的样本，通过几种辐射，原始样本的质量会减少。如果最初有 n 个原子核，经过一段时间后，将有 p 个原子核，其中 p<n。我们知道衰变的量与样本中存在的原子核数量成正比。 −dNdt∝N⟹−dNdt=λN 该符号为负，因为速率随着原子核数量的减少而_减小_。我们将重构此方程，以便我们可以在 Julia 中实现该方程，并使用求解器对其进行求解。 −dNdt∝N⟹−dNdt=λN⟹dNN=−λdt Julia 包 DifferentialEquations.jl 期望函数采用以下格式：

[](https://ritog.github.io/posts/1st-order-DE-julia/<https:/ritog.github.io/posts/1st-order-DE-julia/1st_order_DE_julia.html#cb1-1>)f(du, u, p, t) = p * u

其中 du 是第一个导数，u 是函数，t 是我们想要函数的输出的时间跨度，p 是参数或一组参数。从这里，我们可以将 du 构造为 du[1] = -lambda * u[1]。

导入包

[](https://ritog.github.io/posts/1st-order-DE-julia/<https:/ritog.github.io/posts/1st-order-DE-julia/1st_order_DE_julia.html#cb2-1>)import Pkg
[](https://ritog.github.io/posts/1st-order-DE-julia/<https:/ritog.github.io/posts/1st-order-DE-julia/1st_order_DE_julia.html#cb2-2>)Pkg.add("DifferentialEquations")
[](https://ritog.github.io/posts/1st-order-DE-julia/<https:/ritog.github.io/posts/1st-order-DE-julia/1st_order_DE_julia.html#cb2-3>)Pkg.add("Plots")
[](https://ritog.github.io/posts/1st-order-DE-julia/<https:/ritog.github.io/posts/1st-order-DE-julia/1st_order_DE_julia.html#cb2-4>)
[](https://ritog.github.io/posts/1st-order-DE-julia/<https:/ritog.github.io/posts/1st-order-DE-julia/1st_order_DE_julia.html#cb2-5>)using DifferentialEquations, Plots

Resolving package versions...
No Changes to `~/.julia/environments/v1.11/Project.toml`
No Changes to `~/.julia/environments/v1.11/Manifest.toml`
Resolving package versions...
No Changes to `~/.julia/environments/v1.11/Project.toml`
No Changes to `~/.julia/environments/v1.11/Manifest.toml`

函数定义

[](https://ritog.github.io/posts/1st-order-DE-julia/<https:/ritog.github.io/posts/1st-order-DE-julia/1st_order_DE_julia.html#cb4-1>)function radio_decay!(du, u, p, t)
[](https://ritog.github.io/posts/1st-order-DE-julia/<https:/ritog.github.io/posts/1st-order-DE-julia/1st_order_DE_julia.html#cb4-2>)  lambda = -p
[](https://ritog.github.io/posts/1st-order-DE-julia/<https:/ritog.github.io/posts/1st-order-DE-julia/1st_order_DE_julia.html#cb4-3>)  du[1] = lambda * u[1]
[](https://ritog.github.io/posts/1st-order-DE-julia/<https:/ritog.github.io/posts/1st-order-DE-julia/1st_order_DE_julia.html#cb4-4>)end

radio_decay! (generic function with 1 method)

声明参数

[](https://ritog.github.io/posts/1st-order-DE-julia/<https:/ritog.github.io/posts/1st-order-DE-julia/1st_order_DE_julia.html#cb6-1>)N = [60.0e18]      # 初始条件
[](https://ritog.github.io/posts/1st-order-DE-julia/<https:/ritog.github.io/posts/1st-order-DE-julia/1st_order_DE_julia.html#cb6-2>)
[](https://ritog.github.io/posts/1st-order-DE-julia/<https:/ritog.github.io/posts/1st-order-DE-julia/1st_order_DE_julia.html#cb6-3>)# 参数
[](https://ritog.github.io/posts/1st-order-DE-julia/<https:/ritog.github.io/posts/1st-order-DE-julia/1st_order_DE_julia.html#cb6-4>)lambda = 2.28e6     # 单位 SI 制，s^{-1}
[](https://ritog.github.io/posts/1st-order-DE-julia/<https:/ritog.github.io/posts/1st-order-DE-julia/1st_order_DE_julia.html#cb6-5>)p = lambda
[](https://ritog.github.io/posts/1st-order-DE-julia/<https:/ritog.github.io/posts/1st-order-DE-julia/1st_order_DE_julia.html#cb6-6>)
[](https://ritog.github.io/posts/1st-order-DE-julia/<https:/ritog.github.io/posts/1st-order-DE-julia/1st_order_DE_julia.html#cb6-7>)# 时间跨度
[](https://ritog.github.io/posts/1st-order-DE-julia/<https:/ritog.github.io/posts/1st-order-DE-julia/1st_order_DE_julia.html#cb6-8>)tspan = (0.0, 5e-6)   # 时间跨度非常小，因为半衰期非常小

(0.0, 5.0e-6)

定义 ODE 问题

[](https://ritog.github.io/posts/1st-order-DE-julia/<https:/ritog.github.io/posts/1st-order-DE-julia/1st_order_DE_julia.html#cb8-1>)prob = ODEProblem(radio_decay!, N, tspan, p)

ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
Non-trivial mass matrix: false
timespan: (0.0, 5.0e-6)
u0: 1-element Vector{Float64}:
 6.0e19

求解并绘制解

[](https://ritog.github.io/posts/1st-order-DE-julia/<https:/ritog.github.io/posts/1st-order-DE-julia/1st_order_DE_julia.html#cb9-1>)sol = solve(prob)

retcode: Success
Interpolation: 3rd order Hermite
t: 17-element Vector{Float64}:
 0.0
 4.3859649122807024e-7
 6.41342034530829e-7
 9.724902283236061e-7
 1.2462761785169967e-6
 1.5654877939700708e-6
 1.8726226804630439e-6
 2.196788921416817e-6
 2.5203638024729064e-6
 2.8511887164148543e-6
 3.183703238312926e-6
 3.519781092460714e-6
 3.857529202669002e-6
 4.19723875341719e-6
 4.538177668884617e-6
 4.880273689241893e-6
 5.0e-6
u: 17-element Vector{Vector{Float64}}:
 [6.0e19]
 [2.2085932679488193e19]
 [1.3911043910467977e19]
 [6.538685467576463e18]
 [3.5026563834049787e18]
 [1.6917664046442598e18]
 [8.399157113217204e17]
 [4.011209837867301e17]
 [1.9182281652161018e17]
 [9.022998274997682e16]
 [4.227947963649966e16]
 [1.9650908828375188e16]
 [9.098778059095644e15]
 [4.194139573993325e15]
 [1.927908419090074e15]
 [8.838633947171844e14]
 [6.727182427233061e14]

[](https://ritog.github.io/posts/1st-order-DE-julia/<https:/ritog.github.io/posts/1st-order-DE-julia/1st_order_DE_julia.html#cb11-1>)plot(sol, title="Radioactive Decay", xlabel="Time", ylabel="N")

第二个例子：牛顿冷却定律

正如我们之前在一个例子中看到的，我们知道当物体的温度与周围温度之间的差异相对较高时，物体冷却得更快。当差异较小时，冷却需要_更长_的时间。这个问题是由艾萨克·牛顿本人量化的。让我们来理解它。物体的温度 T 随时间 t 变化，它取决于物体温度与其周围环境（环境温度 Ta）之差。 dTdt∝(T−Ta)⟹dTdt=k(T−Ta) 其中 k 是比例常数，其值取决于材料。现在我们将在 Julia 中实现此 DE，并使用 DifferentialEquations.jl 库对其进行求解。

函数定义

[](https://ritog.github.io/posts/1st-order-DE-julia/<https:/ritog.github.io/posts/1st-order-DE-julia/1st_order_DE_julia.html#cb12-1>)function newtons_cooling!(du, u, p, t)
[](https://ritog.github.io/posts/1st-order-DE-julia/<https:/ritog.github.io/posts/1st-order-DE-julia/1st_order_DE_julia.html#cb12-2>)  k, T_a = p
[](https://ritog.github.io/posts/1st-order-DE-julia/<https:/ritog.github.io/posts/1st-order-DE-julia/1st_order_DE_julia.html#cb12-3>)  du[1] = k * (u[1] - T_a)
[](https://ritog.github.io/posts/1st-order-DE-julia/<https:/ritog.github.io/posts/1st-order-DE-julia/1st_order_DE_julia.html#cb12-4>)end

newtons_cooling! (generic function with 1 method)

期望的格式-

[](https://ritog.github.io/posts/1st-order-DE-julia/<https:/ritog.github.io/posts/1st-order-DE-julia/1st_order_DE_julia.html#cb14-1>)f(du, u, p, t) = p * u

其中 du 是第一个导数，u 是函数，t 是我们想要函数的输出的时间跨度，p 是参数或一组参数。在这里，我们使用就地更新来定义 du，因为它具有内存效率，并且在存在微分方程系统的情况下至关重要。这里我们有两个参数——k 和 T_a，它们之间的关系为 T1=k×(T0−Ta)。这就是我们在函数定义的第 3 行中实现的内容。

声明参数

[](https://ritog.github.io/posts/1st-order-DE-julia/<https:/ritog.github.io/posts/1st-order-DE-julia/1st_order_DE_julia.html#cb15-1>)T = [80.0]        # 物体的温度
[](https://ritog.github.io/posts/1st-order-DE-julia/<https:/ritog.github.io/posts/1st-order-DE-julia/1st_order_DE_julia.html#cb15-2>)k = -0.06         # k，冷却常数，为负，因为是冷却
[](https://ritog.github.io/posts/1st-order-DE-julia/<https:/ritog.github.io/posts/1st-order-DE-julia/1st_order_DE_julia.html#cb15-3>)T_a = 20.0        # 环境温度
[](https://ritog.github.io/posts/1st-order-DE-julia/<https:/ritog.github.io/posts/1st-order-DE-julia/1st_order_DE_julia.html#cb15-4>)tspan = (0.0, 300.0)   # 我们希望追踪函数值的时间段，即物体的温度
[](https://ritog.github.io/posts/1st-order-DE-julia/<https:/ritog.github.io/posts/1st-order-DE-julia/1st_order_DE_julia.html#cb15-5>)p = (k, T_a)       # 定义一个参数元组，符合库的预期

(-0.06, 20.0)

定义 ODE 问题

在这里，我们按照求解器包的要求定义 ODEProblem。

[](https://ritog.github.io/posts/1st-order-DE-julia/<https:/ritog.github.io/posts/1st-order-DE-julia/1st_order_DE_julia.html#cb17-1>)prob = ODEProblem(newtons_cooling!, T, tspan, p)

ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
Non-trivial mass matrix: false
timespan: (0.0, 300.0)
u0: 1-element Vector{Float64}:
 80.0

求解 ODE

[](https://ritog.github.io/posts/1st-order-DE-julia/<https:/ritog.github.io/posts/1st-order-DE-julia/1st_order_DE_julia.html#cb18-1>)sol = solve(prob)

retcode: Success
Interpolation: 3rd order Hermite
t: 19-element Vector{Float64}:
  0.0
  0.18593390178623953
  2.045272919648635
  6.383848011370185
 11.959325821307548
 18.73361986647917
 26.74893491369098
 35.99530039828558
 46.545554755847775
 58.51808266239741
 72.13894261696798
 87.73678833989489
 105.79306069340441
 126.99454945496298
 152.3564523878484
 183.42977158769716
 222.68158609043894
 273.9980369813687
 300.0
u: 19-element Vector{Vector{Float64}}:
 [80.0]
 [79.33435782064105]
 [73.07086984267411]
 [60.9075122329383]
 [49.276503335091164]
 [39.49832238666797]
 [32.054188527309826]
 [26.921535861341653]
 [23.675325765485024]
 [21.792012228502433]
 [20.79153270072668]
 [20.31058542718253]
 [20.10523457704901]
 [20.029607729780903]
 [20.006573373974568]
 [20.0011132329579]
 [20.00018200375497]
 [20.000082850038083]
 [20.00001776796253]

[](https://ritog.github.io/posts/1st-order-DE-julia/<https:/ritog.github.io/posts/1st-order-DE-julia/1st_order_DE_julia.html#cb20-1>)plot(sol,
[](https://ritog.github.io/posts/1st-order-DE-julia/<https:/ritog.github.io/posts/1st-order-DE-julia/1st_order_DE_julia.html#cb20-2>)  xlabel="Time, t",
[](https://ritog.github.io/posts/1st-order-DE-julia/<https:/ritog.github.io/posts/1st-order-DE-julia/1st_order_DE_julia.html#cb20-3>)  ylabel="Temperature, T",
[](https://ritog.github.io/posts/1st-order-DE-julia/<https:/ritog.github.io/posts/1st-order-DE-julia/1st_order_DE_julia.html#cb20-4>)  title="Temperature vs. Time")

我们可以很清楚地看到函数 T 的值如何随时间变化。

结论

我们已经了解了一阶微分方程的非常基础和直观的解释，并且了解了 DifferentialEquations.jl 如何求解 DE，以及如何根据该库期望的格式构造 DE。

致谢

我从大学里一些了不起的老师那里学到了很多关于本教程内容的内容，还有一些了不起的教科书，例如 Riley、Hobson、Bence 的 Mathematical Methods for Physics and Engineering 以及 Arfken、Weber、Harris 的 Mathematical Methods for Physicists。我从创建者 Chris Rackauckas 的许多内容以及这个 YouTube 频道 doggo dot jl 中的一些示例视频中了解了 DE 求解器包。

讨论

如果您阅读了这篇文章，并且发现它有趣或有启发性，请告诉我。我会非常喜欢的。如果您有任何批评、建议或想告诉我任何事情，请添加评论或私下告诉我。在 Fediverse、Hacker News 或 [Twitte