弦月之舞

利用SageMath这个难伺候的主，数值求解高阶微分方程和超越方程。
PS: MathJax渲染数学公式的css加载有点慢，得等等。

R-K法求解边界层的Blasius解

平板边界层方程可以写成这种形式：

$$ \frac{\mathrm d^3 f}{\mathrm d \eta^3} + \frac{1}{2} f \frac{\mathrm d^2 f}{\mathrm d \eta^2} = 0$$

边界条件： $\frac{\mathrm d f}{\mathrm d \eta} \to 1$ as $\eta \to \infty$，$\frac{\mathrm d f}{\mathrm d \eta} = f = 0$ at $\eta = 0$.

由于R-K法只能求解1阶ODE，所以需要换个元，写成方程组的形式：

$$ y_0’ = y_1 $$
$$y_1’ = y_2 $$
$$y_2’’ = -\frac{1}{2} y_0 y_2 $$

参照帮助文档里的ode_solver的教程(https://doc.sagemath.org/html/en/reference/calculus/sage/calculus/ode.html)，调用ode_solver求解，画图即可。不用desolve_system_rk4是因为它的语法写得更奇怪，没看懂，也没调出来。

1
2
3
4
5
6
7
8

T=ode_solver()
eq = lambda t,y:[y[1],y[2],-0.5*y[0]*y[2]]
T.function = eq
T.scale_abs = [1e-4, 1e-4, 1e-5]
T.y_0 = [0,0,0.333]
T.error_rel = 1e-4
T.ode_solve(t_span=[0,10],num_points = 100)
T.plot_solution(i=1)

求超越方程组的数值解

这个例子是毕设中求解André Noth的太阳能飞机设计论文中的总体设计方程的。方程中，翼展b，展弦比AR和总质量m是相互影响的变量，写方程时也会陷入循环引用的圈套。此外，变量名和方程中变量的写法也相当吊诡。最后总之是用隐式解方程的办法解出来了，并批量绘图。

1

find_root(m ==m_tot.subs(AR==10,b==4),0.2,10)

注意find_root和implicit_plot中的变量范围需要多次尝试，不然落在无解的区域就会报错。implicit_plot还不能加图例，到现在也没修好，差评。
点这下载计算脚本。

绘制激波关系

此例来源于空气动力学作业。有些复杂的关系式利用implicit_plot画图未遂，于是用parametric_plot画，但是会缺少解的另一半支，这使得最后的结果并不优美。
绘图脚本

至于图长什么样，还是留待SageMath里跑一下吧:)。
(虽然这玩意只有在Linux下才比较好装。)简单的程序也可以用Sage Cell Server跑，至少解ODE的例子是可以跑的。