微分方程式を行列で差分化して解く　その３

今回も前回に引き続き行列を使用して1階微分の微分方程式を解いてみます。前回までは初期値が $$ y(0) = 0 $$ でしたが、そうではない場合について考えていきたいと思います。

前回までの内容は以下になります。

otepipi.hatenablog.com

今回解く数式は以下になります。

$$ \begin{equation} \frac{dy}{dx}=3x^{2}, y(0) = 10 \end{equation} $$

この方程式の解析解は $$ y = x^{3} + 10 $$ になります。

この場合も前回同様のやり方（⊿x = 0.01)で解けるかどうか確認してみます。

数式としては以下のようになります。

$$ \begin{bmatrix} 1 & 0 & \cdots & 0 & 0 \\ -1 & 1 & \cdots & 0 & 0 \\ 0 & -1 & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \ddots & \vdots \\ 0 & 0 & \cdots & -1 & 1 \end{bmatrix} \begin{bmatrix} y(0) \\ y(0.01) \\ y(0.02) \\ \vdots \\ y(3) \end{bmatrix} 　= 0.01 \times \begin{bmatrix} 10\\ 0.0003 \\ 0.0012 \\ \vdots \\ 27 \end{bmatrix} $$

右辺の1行目を10に変更しています。

ではこの解法をJuliaで実装してみます。

using LinearAlgebra
using Gadfly

dv = ones(301)
ev = -1 *ones(300)
A_new =  Bidiagonal(dv,ev, :L)

k = 0:0.01:3
b_new = reshape(0.01*3*k.^2,301,1)

b_new[1] = 10

y_new = A_new \ b_new

println(y_new)

x = 0:0.01:3
y_real = x.^3 .+ 10
plot(
    layer(x = x, y = y_real,Geom.line,Theme(default_color=colorant"red")),
    layer(x = x, y = y_new)
    )