在 Julia 中实现 Runge-Kutta 方法

我在 Julia 中实现了 Runge-Kutta 方法。

在许多编程语言中，函数也是函数的参数，可以在函数中直接引用。
feval(func, t(k), Y(:,k)) 可以简单地写成func(t(k), Y(:,k))。

## ルンゲ・クッタ法
function ClassicalRungeKutta(func, t, Y0)
  # func:微分方程式の関数ハンドル
  # t：積分範囲（hで分割済み）
  # Y0：初期値
  nrow = length(Y0)
  ncol = length(t)
  Y = zeros(nrow, ncol)
  Y[:, 1] = Y0 # 初期値の代入
  h = t[2] - t[1] # 刻み幅
  for k = 1:ncol-1
      k1 = func(t[k], Y[:, k])
      k2 = func(t[k] + h/2, Y[:, k] + h*k1/2)
      k3 = func(t[k] + h/2, Y[:, k] + h*k2/2)
      k4 = func(t[k] +h, Y[:, k] + h*k3)
      Y[:, k+1] = Y[:, k] + h*(k1 + 2*k2 + 2*k3 + k4)/6
  end
  Y
end

## 斜方投射（鉛直上向き）の運動方程式
function MotionEquation(t, Y)
    g = 9.8
    nrow = length(Y)
    Ydot = zeros(nrow, 1)
    Ydot[1, 1] = Y[2]
    Ydot[2, 1] = -g
    Ydot
end

# 質点の斜方投射

##　物理パラメータと運動方程式の初期条件

x0 = 0; y0 = 0 # 原点
g = 9.8 # 重力加速度
theta = (45/180)*pi # 45度の方向へ投げる
v0 = 10 # 初速度 m/s

# 運動方程式の初期条件
Y0 = [y0 v0*sin(theta)]

# シミュレーション時間
h = 10^-4 #刻み幅
t = 0:h:5

## 　運動方程式を解いてボールの軌跡を描画

x = @. x0 + v0*cos(theta)*t;
y = ClassicalRungeKutta(MotionEquation, t, Y0) # 近似解
yTrue = @. y0 + v0*sin(theta)*t - (g*t.^2)/2; # 真の解

using Plots
gr(fontfamily="IPAMincho", foreground_color_legend = nothing)
plot(x, y[1, :], xlims=(0, 10.1), ylims=(0, 3), label="真の解")
plot!(x, yTrue, label="近似解", title="斜方投射", xlabel="x[m]", ylabel="y[m]")

Julia でルンゲ・クッタ法を実装してみる

原创声明：本文系作者授权爱码网发表，未经许可，不得转载;

原文地址：https://www.likecs.com/show-308628559.html