Julia 求解微分方程

发表于 分类于 ， 阅读次数： 本文字数： 924 阅读时长 ≈ 2 分钟

本文为使用 DifferentialEquations.jl 求解微分方程的一个例子。

安装求解器

1
]add DifferentialEquations

为方便绘图，建议安装 Plots.jl 这个包

一个例子

在这个实例中，我们将求解方程
$$
\begin{cases}
\frac{\rm du}{\rm dt}=1.01u\\
u(0)=0.5
\end{cases}
$$

定义问题

1
2
3
4
5
6
7
function eq!(du, u, p, t)
    du = p*u
end
u0 = 0.5
p = 1.01  //传入的参数
tspan = (0.0,1.0)  //求解的时间区间
prob = ODEProblem(eq!,u0,tspan,p)

求解及求解器参数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
//直接求解
sol = solve(prob)

//相对误差控制在1e-6以内
sol = solve(prob,reltol=1e-6)

//每隔0.1个时间点保存一次，此参数 也可传入一个要保存点的数组
sol = solve(prob,saveat=0.1)

//只保存端点值
sol = solve(prob,save_everystep=false)

//必须求解1.0, 1.5, 2.0这几个点，此参数可以与 saveat 连用，也可以帮助处理方程不连续
sol = solve(prob,tstops=[1.0, 1.5, 2.0])

//插值只使用一阶插值，默认为三阶，因为有一些阶三阶导数不连续
sol = solve(prob,dense=false)

绘图

绘图需要使用 Plots.jl 这个包

1
2
3
4
5
6
7
8
9
10
11
12
13
14
//直接绘图
plot(sol)

//只使用一阶插值绘图，默认为三阶插值
plot(sol,denseplot=false)

//只绘制0-1s这个时间段
plot(sol, tspan=(0.0,1.0))

//只绘制解中第一个变量，这对于微分方程组有用
plot(sol,vars=[1])

//绘制之后保存图片
savefig("myplot.png")

查看解

1
2
//获得在t处的解值
sol(t)

更多例子及用法

参考官方文档：https://diffeq.sciml.ai/stable/