0

変数を式に代入する方法について。

D=0で計算していますが、D=[0,0.01,0.1]の三つをそれぞれf(x,v,t)に代入し、1つのグラフで比較するにはどうすればいいですか?

import numpy as np
import matplotlib.pyplot as plt

def f(x, v, t):
    M = 1.0
    E = 0.38 if (t <= 0.15) else 1.2
    J = (2*0.8)/(120*np.pi)
    D = 0
    return (M-E*np.sin(x)-D*v)/J

t0 = 0.0
t1 = 5.0

N = 500

del_t = (t1-t0)/N # time grid

tpoints = np.arange(t0, t1, del_t)
xpoints = []
vpoints = []

# initial condition
x0 = 0.653
v0 = 0.0

x, v = x0, v0
for t in tpoints:
    xpoints.append(x)
    vpoints.append(v)
    k1v =f(x, v, t)*del_t
    k1x = v * del_t

    k2v =  f(x+k1x/2, v+k1v/2, t+del_t/2 )*del_t
    k2x =(v+k1v/2)*del_t 

    k3v =f (x+k2x/2, v+k2v/2, t+del_t/2 )*del_t
    k3x =(v+k2v/2 ) *del_t 

    k4v = f(x+k3x, v+k3v, t+del_t )*del_t
    k4x = (v+k3v )*del_t 

    v += (k1v + 2 * k2v + 2* k3v + k4v)/6
    x += (k1x + 2 * k2x + 2* k3x + k4x)/6

plt.plot (tpoints, xpoints, 'o',label='4th order Runge-Kutta')
plt.xlabel("t",  fontsize=24)
plt.ylabel("x(t)",  fontsize=24)

plt.show()

2 件の回答 2

0

【修正方針】
今あるコードを最大限に活かした変更とします。

【修正方法】
まず、

t0 = 0.0

から

x += (k1x + 2 * k2x + 2* k3x + k4x)/6

までを関数化します。

変数Dtpointsxpointsをglobal変数とします。※関数fD=0は削除

後は、Dを変えながら関数化した処理を呼び出し、グラフを表示するだけです。

【修正後のコード】

import numpy as np
import matplotlib.pyplot as plt

def f(x, v, t):
    global D
    M = 1.0
    E = 0.38 if (t <= 0.15) else 1.2
    J = (2*0.8)/(120*np.pi)
    return (M-E*np.sin(x)-D*v)/J

def func():
    global tpoints
    global xpoints
    t0 = 0.0
    t1 = 5.0

    N = 500

    del_t = (t1-t0)/N # time grid

    tpoints = np.arange(t0, t1, del_t)
    xpoints = []
    vpoints = []

    # initial condition
    x0 = 0.653
    v0 = 0.0

    x, v = x0, v0
    for t in tpoints:
        xpoints.append(x)
        vpoints.append(v)
        k1v =f(x, v, t)*del_t
        k1x = v * del_t

        k2v =  f(x+k1x/2, v+k1v/2, t+del_t/2 )*del_t
        k2x =(v+k1v/2)*del_t 

        k3v =f (x+k2x/2, v+k2v/2, t+del_t/2 )*del_t
        k3x =(v+k2v/2 ) *del_t 

        k4v = f(x+k3x, v+k3v, t+del_t )*del_t
        k4x = (v+k3v )*del_t 

        v += (k1v + 2 * k2v + 2* k3v + k4v)/6
        x += (k1x + 2 * k2x + 2* k3x + k4x)/6

Ds = [0, 0.01, 0.1]
for d in Ds:
    D = d
    func()
    plt.plot (tpoints, xpoints, 'o',label='4th order Runge-Kutta')
    plt.xlabel("t",  fontsize=24)
    plt.ylabel("x(t)",  fontsize=24)

plt.show()

【実行結果】

【注意】
 グラフのスケールを調整しないと、Dが0.1のときの微妙な曲線がなだらかになってしまいます。

0

関数 f のパラメーターに d を追加して、Runge-Kutta 法で計算処理を行っている部分を関数化します(rk4)。

import numpy as np
import matplotlib.pyplot as plt

M = 1.0
J = (2*0.8)/(120*np.pi)

def f(x, v, t, d):
    E = 0.38 if (t <= 0.15) else 1.2
    return (M-E*np.sin(x)-d*v)/J

def rk4(x0, v0, tpoints, d, delta):
  x, v = x0, v0
  xpoints, vpoints = [], []
  for t in tpoints:
    xpoints.append(x)
    vpoints.append(v)

    k1v = f(x, v, t, d)*delta
    k1x = v*delta

    k2v = f(x+k1x/2, v+k1v/2, t+delta/2, d)*delta
    k2x = (v+k1v/2)*delta 

    k3v = f(x+k2x/2, v+k2v/2, t+delta/2, d)*delta
    k3x = (v+k2v/2)*delta 

    k4v = f(x+k3x, v+k3v, t+delta, d)*delta
    k4x = (v+k3v)*delta 

    x += (k1x + 2 * k2x + 2* k3x + k4x)/6
    v += (k1v + 2 * k2v + 2* k3v + k4v)/6

  return xpoints, vpoints

if __name__ == '__main__':

  t0, t1, N = 0.0, 5.0, 500
  del_t = (t1-t0)/N # time grid
  tpoints = np.arange(t0, t1, del_t)

  # initial condition
  x0, v0 = 0.653, 0.0
  D = (0, 0.01, 0.1)
  xvpoints = [rk4(x0, v0, tpoints, d, del_t) for d in D]
  vlabels = [f'D = {d}' for d in D]

  for xv, label in zip(xvpoints, vlabels):
    plt.plot(tpoints, xv[0], 'o', label=label)
  plt.legend(loc='upper left', fontsize=10)
  plt.title('4th order Runge-Kutta')
  plt.xlabel("t",  fontsize=24)
  plt.ylabel("x(t)",  fontsize=24)

  plt.show()

スケールの差異が大きすぎるので全体的な傾向の違いしか判りませんけれども。

この質問に回答するには、ログインする必要があります。

求めていた回答ではありませんか? のタグが付いた他の質問を参照する。