ModelingToolkit建模方法

符号计算内涵

符号计算的内涵很简单，其本质就是处理、简化方程的未知数。例如: $(x+y)z = xy + yz$ 又例如: $\frac{dy}{dt} = 2t => y = t^2 + C $ 它们本质上都是数学符号的推演，而没有涉及到数字的计算。

一般的情况下，可以计算问题都是人为地把计算推演到数字计算的层面，中间的符号推演、化简的过程都不由计算机处理。随着计算机能力的日益强大，以及计算理论的发展。计算机也能开始进行符号计算。

例如，$(x+y)z = xy + yz$就完全能由计算机计算并得到结果。

MTK，就是这样一种符号计算包。在DE中，我们构建的是方程（function），通过方程来建立模型。而在MTK中，我们构建的是符号，通过符号建立模型。

MTK符号计算实例

考虑洛伦兹方程：

\[\frac{dx}{dt} = \sigma(y-x) \\ \frac{dy}{dt} = x(\rho-z)-y \\ \frac{dz}{dt} = xy - \beta z\]

通过MTK构建：

using ModelingToolkit
using DifferentialEquations

@variables t u(t)[1:3]
@parameters p[1:3]

D = Differential(t)
eqs =[
    D(u[1]) ~ p[1]*(u[2]-u[1])
    D(u[2]) ~ u[1]*(p[2]-u[3]) - u[2]
    D(u[3]) ~ u[1]*u[2] - p[3]*u[3]
]
@named Sys = ODESystem(eqs)

\[ \begin{align} \frac{\mathrm{d}}{\mathrm{d}t} u(t)_1 =& \left( - u(t)_1 + u(t)_2 \right) p_1 \\ \frac{\mathrm{d}}{\mathrm{d}t} u(t)_2 =& - u(t)_2 + \left( - u(t)_3 + p_2 \right) u(t)_1 \\ \frac{\mathrm{d}}{\mathrm{d}t} u(t)_3 =& u(t)_1 u(t)_2 - p_3 u(t)_3 \end{align} \]

其中：

• t是时间变量
• u是关于时间的依赖变量向量
• 参数向量p
• 微分符号D
• eqs本质是一个描述符号方程的数组

@variables定义了符号变量，@parameters p[1:3]定义了参数。方程由ODESystem结构体储存，可以看到其中描述“=”的符号为“~”。这些表明，上述代码在处理符号，这是一套符号描述体系！

因为是处理符号时，当然也可以更加具象，将符号定义为$x,y,z$

using ModelingToolkit
using DifferentialEquations

@variables t x(t) y(t) z(t)
@parameters  σ ρ β
D = Differential(t)
eqs =[
    D(x) ~ σ*(y-x)
    D(y) ~ x*(ρ-z) - y
    D(z) ~ x*y - β*z
]
@named sys = ODESystem(eqs)

\[ \begin{align} \frac{\mathrm{d} x\left( t \right)}{\mathrm{d}t} =& \sigma \left( - x\left( t \right) + y\left( t \right) \right) \\ \frac{\mathrm{d} y\left( t \right)}{\mathrm{d}t} =& - y\left( t \right) + \left( \rho - z\left( t \right) \right) x\left( t \right) \\ \frac{\mathrm{d} z\left( t \right)}{\mathrm{d}t} =& x\left( t \right) y\left( t \right) - \beta z\left( t \right) \end{align} \]

由于符号运算的展现形式更加贴近自然语言，我们几乎可以不加处理地建立与原数学方程几乎一样的符号方程描述，这是符号计算的优势之一。

在明确MTK的系统构建之后，与DifferentialEquations一样，最后需要明确定义的问题以及初值等要素。

tspan = (0.0,100.0)
u0 =[
    x => 1.0
    y => 0.0
    z => 0.0
]
p=[
    σ => 10.0
    ρ => 28.0
    β => 8/3
]
prob = ODEProblem(sys,u0,tspan,p)
sol = solve(prob,Tsit5())

retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 1263-element Vector{Float64}:
   0.0
   3.5678604836301404e-5
   0.0003924646531993154
   0.0032624077544510573
   0.009058075635317072
   0.01695646895607931
   0.02768995855685593
   0.04185635042021763
   0.06024041165841079
   0.08368541255159562
   ⋮
  99.30760258626904
  99.39665422328268
  99.49536147459878
  99.58822928767293
  99.68983993598462
  99.77864535713971
  99.85744078539504
  99.93773320913628
 100.0
u: 1263-element Vector{Vector{Float64}}:
 [1.0, 0.0, 0.0]
 [0.9996434557625105, 0.0009988049817849058, 1.781434788799208e-8]
 [0.9961045497425811, 0.010965399721242457, 2.146955365838907e-6]
 [0.9693591634199452, 0.08977060667778931, 0.0001438018342266937]
 [0.9242043615038835, 0.24228912482984957, 0.0010461623302512404]
 [0.8800455868998046, 0.43873645009348244, 0.0034242593451028745]
 [0.8483309847495312, 0.6915629321083602, 0.008487624590227805]
 [0.8495036669651213, 1.0145426355349096, 0.01821208962127994]
 [0.9139069574560097, 1.4425599806525806, 0.03669382197085303]
 [1.088863826836895, 2.052326595543049, 0.0740257368585531]
 ⋮
 [4.669609096878053, 3.061564434452441, 25.1424735017959]
 [4.188801916573263, 4.617474401440693, 21.09864175382292]
 [5.559603854699961, 7.905631612648314, 18.79323210016923]
 [8.556629716266505, 12.533041060088328, 20.6623639692711]
 [12.280585075547771, 14.505154761545633, 29.332088452699942]
 [11.736883151600804, 8.279294641640229, 34.68007510231878]
 [8.10973327066804, 3.2495066495235854, 31.97052076740117]
 [4.958629886040755, 2.194919965065022, 26.948439650907677]
 [3.8020065515435855, 2.787021797920187, 23.420567509786622]

与DE不同是，MTK传值时需要指定具体变量（参数的值）。 最后能得到与DE计算同样的结果。

全部代码：

using ModelingToolkit
using DifferentialEquations
@variables t x(t) y(t) z(t)
@parameters  σ ρ β
D = Differential(t)
eqs =[
    D(x) ~ σ*(y-x)
    D(y) ~ x*(ρ-z) - y
    D(z) ~ x*y - β*z
]
@named sys = ODESystem(eqs)
tspan = (0.0,100.0)
u0 =[
    x => 1.0
    y => 0.0
    z => 0.0
]
p=[
    σ => 10.0
    ρ => 28.0
    β => 8/3
]
prob = ODEProblem(sys,u0,tspan,p)
sol = solve(prob,Tsit5())
using Plots
plot(sol, idxs=(x,y,z))

DE与MTK的对比

同一问题，可用MTK与DE来求解。但它们并不完全等价，异同是客观存在的：

• 不同之处： ModelingToolkit基于符号体系描述问题，等价于在DifferentialEquations问题描述的上层添加了一层符号系统，通过求解符号系统之后，能得到真正的微分方程问题。

• 相同之处： 它们最终得到的问题形式是一样的，比如说最终都会生成ODEProblem，也就是说它们可以通过同一个求解器求解。

MTK凭借着符号计算，在模型建立时会更加方便，更加友好。符号推演功能相当于减轻了模型构造时的工作量，即在数学模型到计算机可计算的模型这个过程中，计算机又往前前进了一步。所以，MTK从某种程度上来说会比DE好用很多。

从另一个角度来看，DE是MTK的底层求解器，意味着其比MTK更加直接，速度更快。MTK在DE基础之上构建了符号计算系统，其形式与功能上让人惊叹。但要明白的是，这种提供便捷的符号计算系统是要使用资源的，MTK需要花费更多资源去计算、化简符号。好用意味着大量的资源消耗在用户友好方面，效率高意味着资源消耗在真正的问题解决上。这是软件发展过程中始终存在的矛盾。

关于DE与MTK的取舍，根据问题的特性选择即可。

此外，MTK的符号系统还有更大的妙用，在后面的章节中会逐渐展开。