Julia 求解ODE方程 2

为什么要使用Julia？ (这些都是没用的话）

很多人觉得语言只是工具，没有本质区别。其实并不是这样，举个例子，你想编一个大型的程序，其中你只关心你提出的一个方法，而其他大部分功能和方法都只是为了实现，并不是你的创新点，你也没必要真的去学习它们，只需要了解到它们各有什么优点和在做什么就可以了。这种时候就体现出了高级语言站在巨人的肩膀上的优势了，有极多的库和包可以用，甚至你的程序百分九十的功能都是基于别人封装好的包来做的。 这有什么优点呢，第一，节约了你大部分时间，没必要编写和调试本来无需去做的工作，第二，成熟的包和库都经过了大量人力的优化，计算速度和稳定性远远强于你自己短期编的程序，第三，如果是开源的库，如果你有更多的需要，可以在原框架上做自己的修改。

当然，Fortran 和c等初级语言也有现成的代码可以用，比如MKL等商用包，但是用起来真的超级不方便，比如你想实现一个矩阵求逆的运算，需要下载、安装等等一系列繁琐的工作，然后调用时候也很不方便，比如会同时需要给十几个甚至几十个参数输入。另外调用了别人库的程序可移植性变得很差，经常换个系统，换个电脑就无法调试通。 所以大家仍然喜欢自己去编绝大部分功能。 当然，如果是在学习阶段，充分理解某个计算流程，初级语言仍然是很有必要的。

这是高级语言，包括julia的一个优点。 但是，为什么我没有选择Python或者matlab、octave等呢，那就是计算速度问题了， 我将要做的工作大部分和大型计算相关，目前代码基本上都是Fortran和c的，在计算速度方面，FOrtran基本上可以秒空气了，然而，Julia也是为了科学计算开发的，计算速度并不弱于Fortran和c太多，另外，由于Fortran编写中，如果你不考虑优化的话，甚至还没有julia程序快。即便牺牲一点点计算速度，换取更多的个人时间都是非常值得的。

另外julia无需胶水语言就可以编译Fortran，C，Python等代码，代码的继承能力也是极强的。不断完善的各类高性能计算、机器学习的库，以及并行和Gpu等极多的优点，还有更高的可读性，释义型语言的优势等等，都非常具有吸引力。。。

当然，julia上手难度还是颇高的，尤其对于没有任何高级语言经验，只会面向过程编程，需要学习和练习很多，但是我觉得这些都是值得的。另外，julia毕竟只是个才发展10年的语言，目前还缺乏足够多的应用，尤其国内几乎没有公司，只有极少数的个人在使用julia，仍然算是非常小众的东东，但是作为个人应用，真的足够了。 非常看好它的前景。

好了言归正传，在给一个DifferentialEquations求解ODEs方程组的例子，

---

Lorenz 方程

$$\begin{aligned} &\frac{dx}{dt} = \sigma (y-x)\\ &\frac{dy}{dt}=x(\rho -z)-y\\ &\frac{dz}{dy}=xy-\beta z \end{aligned}$$
同样我们需要先定义函数

    function larenz(du,u,p,t)
        du[1]= 10.0(u[2]-u[1])
        du[2] = u[1] * (28.0-u[3]) -u[2]
        du[3]= u[1]*u[2] - (8/3)u[3]
    end

当然这么写很不好，换一套系数还得改程序不是？

    function larenz(du,u,p,t)
        du[1]= p[1]*(u[2]-u[1])
        du[2] = u[1] * (p[2]-u[3]) -u[2]
        du[3]= u[1]*u[2] - p[3]*u[3]
    end

这样就好多了，然而DE包还给了一个宏可以用，很骚的操作。

        g = @ode_def LorenzExample begin
            dx = σ*(y-z)
            dy = x*(ρ -z)- y
            dz = x*y -β*z
        end σ ρ β

是不是骚气多了，这还不算完，调用ParameterizeFunctions 包还可以实现更骚的操作

g = @odedef Lorenz2 begin
    df₁ = σ*(f₂ - f₃)
    df₂ = f₁*(ρ - f₃) - f₂
    df₃ = f₁*f₂ - β*z
end σ ρ β

如果符号已经不能表达你绝望的写代码的心情了

那我们就用emoji吧。当然这有点太浮夸了，我们还是采用g来计算，
第二步，给出初始条件和时间、参数等

g1 = @ode_def Lorenz2 begin
    df₁ = σ*(f₂ - f₃)
    df₂ = f₁*(ρ - f₃) - f₂
    df₃ = f₁*f₂ - β*f₃
end σ ρ β

tspan=(0.0,1.0)
p=[10.0,28.0,8/3]
u0=[1.0,0.0,0.0]

prob=ODEProblem(g1,u0,tspan,p)
sol = solve(prob)
using Plots
plot(sol,vars=(1,2,3))

最后画个三维图，

当然这里Plots对字符支持还是有点问题，有时间我再研究下画图的包吧，这个也是非常有必要的，太晚了，睡觉zzz 明天有时间再说。