66865


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

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

[code]

    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()
    
[/code]


CC BY-SA 4.0
2020-05-22T01:09:44.307
66855
2020-05-22T07:31:43.380
2020-05-22T01:51:31.590
19110
40181
question
python matplotlib 数値解析
グラフを複数同時に出したい


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

**【修正方法】**  
まず、

> t0 = 0.0

から

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

までを関数化します。

変数`D`、`tpoints`、`xpoints`をglobal変数とします。※関数`f`の`D=0`は削除

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

**【修正後のコード】**

[code]

    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()
    
[/code]

**【実行結果】**  
[![](https://i.stack.imgur.com/ooBWM.png)](https://i.stack.imgur.com/ooBWM.png)

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


CC BY-SA 4.0
2020-05-22T03:15:00.513
66865
2020-05-22T07:31:43.380
2020-05-22T07:31:43.380
35558
35558
66855
answer


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

[code]

    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()
    
[/code]

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

[![](https://i.stack.imgur.com/BZkUk.png)](https://i.stack.imgur.com/BZkUk.png)


CC BY-SA 4.0
2020-05-22T03:40:58.943
66866
2020-05-22T03:40:58.943
66855
answer