# ode89
求解非刚性微分方程 - 高阶方法
函数库: TyMath
# 语法
t,y, = ode89(odefun,tspan,y0,options)
t,y,te,ye,ie = ode89(odefun,tspan,y0,options)
sol = ode89(___;output_sol=true)
# 说明
t,y, = ode89(odefun,tspan,y0)（其中 tspan = [t0 tf]）求微分方程组
所有 ODE 求解器都可以解算 y′=f(t,y) 形式的方程组，或涉及质量矩阵 M(t,y)y′=f(t,y) 的问题。求解器都使用类似的语法。ode23s 求解器只能解算质量矩阵为常量的问题。ode15s 和 ode23t 可以解算具有奇异质量矩阵的问题，称为微分代数方程 (DAE)。使用 odeset 的 :mass 选项指定质量矩阵。 示例
t,y, = ode89(odefun,tspan,y0,options) 还使用由 options（使用 odeset 函数创建的参数）定义的积分设置。例如，使用 :abstol 和 :reltol 选项指定绝对误差容限和相对误差容限，或者使用 :mass 选项提供质量矩阵。
t,y,te,ye,ie = ode89(odefun,tspan,y0,options) 还求 (t,y) 的函数（称为事件函数）在何处为零。在输出中，te 是事件的时间，ye 是事件发生时的解，ie 是触发的事件的索引。
对于每个事件函数，应指定积分是否在零点处终止以及过零方向是否重要。为此，请将 :events 属性设置为函数（例如 myEventFcn），并创建一个对应的函数： myEventFcn = (t,y) -> (value,isterminal,direction)。
sol = ode89(___;output_sol=true) 返回一个结构体，您可以将该结构体与 deval 结合使用来计算区间 [t0 tf] 中任意点位置的解。您可以使用上述语法中的任何输入参数组合。
# 示例
具有一个解分量的 ODE
在对求解器的调用中，可将只有一个解分量的简单 ODE 指定为匿名函数。该匿名函数必须同时接受两个输入 (t,y)，即使在该函数中一个输入未使用也是如此。
解算 ODE
指定时间区间 [0 5] 和初始条件 y0 = 0。
using TyMath
using TyPlot
tspan = [0 5];
y0 = 0;
t,y, = ode89((t,y)->2*t, tspan, y0);
对解绘图
plot(t,y,"-o")
解算非刚性方程
van der Pol 方程为二阶 ODE
其中
函数 vdp1 代表使用
using TyMath
using TyPlot
function vdp1(t,y)
dydt = [y[2]; (1-y[1]^2)*y[2]-y[1]]
return dydt
end
使用 ode89 函数、时间区间 [0 20] 和初始值 [2 0] 来解算该 ODE。生成的输出即为时间点 t 的列向量和解数组 y。y 中的每一行都与 t 的相应行中返回的时间相对应。y 的第一列与
t,y, = ode89(vdp1,[0 20],[2; 0]);
绘制
plot(t,y[:,1],"-o",t,y[:,2],"-o")
title(raw"$Solution of van der Pol Equation (\mu = 1) with ode89$");
xlabel("Time t");
ylabel("Solution y");
legend(raw"$y_1$",raw"$y_2$")
向 ODE 函数传递额外的参数
ode89 仅适用于使用两个输入参数（t 和 y）的函数。但是，通过在函数外部定义参数并在指定函数句柄时传递这些参数，可以传入额外参数。
解算 ODE
将该方程重写为一阶方程组可以得到
此示例部函数 odefcn 将此方程组表示为接受四个输入参数（t、y、A 和 B）的函数。
using TyMath
using TyPlot
function odefcn(t, y, A, B)
dydt = [y[2]; (A/B)*t.*y[1]]
return dydt
end
使用 ode89 解算 ODE。指定函数句柄，使其将 A 和 B 的预定义值传递给 odefcn。
A = 1;
B = 2;
tspan = [0 5];
y0 = [0 0.01];
t,y, = ode89((t,y)->odefcn(t,y,A,B), tspan, y0);
绘制结果。
plot(t,y[:,1],"-o",t,y[:,2],"-.")
具有严格误差容限的 ODE
与 ode45 相比，ode113、ode78 和 ode89 求解器可以更好地求解具有严格误差容限的问题。这些求解器更胜一筹的常见情况是在轨道动力学问题中，这些问题的解曲线是平滑的且在求解器的每个时间步中都需要很高的准确度。
二体问题涉及两个相互作用的质点 m1 和 m2 在一个公共平面上沿轨道运动。在本示例中，一个质点明显大于另一个。以较重的天体为原点，运动方程为：
其中
要求解该问题，首先使用以下替换公式转换为包含四个一阶 ODE 的方程组：
替换之后生成的一阶方程组为：
函数 twobodyode 为二体问题的方程组进行编码。
using TyMath
using TyPlot
function twobodyode(t, y)
r = sqrt(y[1]^2 + y[3]^2)
dy = [y[2]; -y[1] / r^3; y[4]; -y[3] / r^3]
return dy
end
使用 ode89 求解 ODE 方程组。对 RelTol 指定严格误差容限 1e-13，对 AbsTol 使用严格误差容限 1e-14。
opts = odeset(; reltol=1e-13, abstol=1e-14, stats=:on)
tspan = [0 10 * pi]
y0 = [2 0 0 0.5]
t, y, = ode89(twobodyode, tspan, y0, opts);
243 个成功步骤
0 次失败尝试
5103 次函数计算
plot(t,y)
legend(["x","x'","y","y'"],loc="southeast")
title("Position and Velocity Components")
figure()
plot(y[:,1],y[:,3],"-o",0,0,"ro")
axis("equal")
title("Orbit of Smaller Mass")
与 ode45 相比，ode89 求解器能够通过更少的步骤和函数计算更快地获得解。
# 输入参数
odefun - 要求解的函数函数
要求解的函数，指定为指向待积分器的句柄。
对于标量 t 和列向量 y 来说，函数 odefun = (t,y) -> dydt 必须返回列向量 dydt，该列向量对应于 f(t,y)。odefun 必须同时接受输入参数 t 和 y，即使其中一个参数未在函数中使用也是如此。
例如，要解算
function odefun(t,y)
return 5*y .- 3
end
对于方程组，odefun 的输出为向量。向量中的每个元素是一个方程的解。例如，要求解
使用函数：
function odefun(t,y)
dydt = [y[1]+2*y[2];3*y[1]+2*y[2]]
return dydt
end
示例： myFcn
数据类型： Function
tspan - 积分区间向量
积分区间，指定为向量。至少，tspan 必须是指定初始时间和最终时间的二元素向量 [t0 tf]。要获取 t0 到 tf 之间的特定时间的解，请使用 [t0,t1,t2,...,tf] 形式的长向量。tspan 中的元素必须单调递增或单调递减。
求解器在初始时间 tspan[1] 施加由 y0 给出的初始条件，然后求 tspan[1] 到 tspan[end] 的积分：
如果 tspan 有两个元素 [t0 tf]，求解器将返回在该区间内的每个内部积分步计算的解；
如果 tspan 包含两个以上的元素，[t0,t1,t2,...,tf]，求解器将返回在给定点处计算的解。但是，求解器不会精确步进到 tspan 中指定的每个点。此时，求解器使用自己的内部积分步来计算解，然后在 tspan 中请求的各点处计算解。在指定点处生成的解与在每个内部积分步计算的解具有相同的准确度级别。
指定多个中间点对计算效率影响甚微，但对于大型系统可能会影响内存管理。
求解器使用 tspan 的值计算 InitialStep 和 MaxStep 的合适值：
如果 tspan 包含多个中间点 [t0,t1,t2,...,tf]，则指定的点表示了问题的规模，这可能影响求解器使用的 :initialstep 的值。因此，根据您是将 tspan 指定为二元素向量还是包含中间点的向量，求解器获得的解可能有所不同；
tspan 中的初始值和最终值用于计算最大步长 :maxstep。因此，更改 tspan 中的初始值或最终值可能导致求解器使用不同步长序列，从而可能会更改解。
示例： [1 10]
示例： [1 3 5 7 9 10]
数据类型： Int | Float
y0 - 初始条件向量
初始条件，指定为向量。y0 的长度必须与 odefun 的向量输出相同，使 y0 为 odefun 中定义的每个方程包含一个初始条件。
数据类型： Int | Float
options - options 结构体OdeOption
options 结构体，指定为结构体数组。使用 odeset 函数创建或修改 options 结构体。
示例： options = odeset(reltol=1e-5,stats=:on) 指定 1e-5 的相对误差容限、打开求解器统计信息的显示。
数据类型： OdeOption
# 输出参数
t - 求值点列向量
计算点，以列向量形式返回。
如果 tspan 包含两个元素 [t0 tf]，则 t 包含用于执行积分的内部计算点；
如果 tspan 包含两个以上元素，则 t 与 tspan 相同。
数据类型： Int | Float
y - 解数组
解，以数组形式返回。y 中的每一行都与 t 的相应行中返回的值处的解相对应。
数据类型： Int | Float
te - 事件的时间列向量
事件的时间，以列向量形式返回。te 中的事件时间对应于 ye 中返回的解，而 ie 指定发生了哪个事件。
数据类型： Int | Float
ye - 事件时间处的解数组
事件时间的解，以数组形式返回。te 中的事件时间对应于 ye 中返回的解，而 ie 指定发生了哪个事件。
数据类型： Int | Float
ie - 触发的事件函数的索引列向量
触发的事件函数的索引，以列向量形式返回。te 中的事件时间对应于 ye 中返回的解，而 ie 指定发生了哪个事件。
数据类型： Int
sol - 用于计算的结构体结构体
用于计算的结构体，以结构体数组形式返回。
| 结构体字段 | 说明 |
|---|---|
| sol.x | 求解器选择的步的行向量。 |
| sol.y | 解。每列 sol.y[:,i] 包含时间 sol.x[i] 处的解。 |
| sol.solver | 求解器名称。 |
此外，如果指定了 odeset 的 :events 选项并且检测到事件，则 sol 的下列字段不为nothing：
| 结构体字段 | 说明 |
|---|---|
| sol.xe | 事件发生的点。sol.xe[end] 包含终止事件（如果有）的确切点。 |
| sol.ye | 与 sol.xe 中的事件相对应的解。 |
| sol.ie | :events 选项中指定的函数所返回的向量的索引。这些值指示求解器检测到的事件。 |
# 算法
ode89 是 Verner 的“最稳健”的 8 阶连续外推 Runge-Kutta 9(8) 法的一种实现。解的精度提升到了 9 阶。8 阶连续外推需要对 odefun 进行五次额外计算，但仅限于需要插值的步。