微分方程式を解きたいです．

Vrを0.5から24まで変化させ，各Vrでもｔを0から8秒間（サンプリング点数：16384点）の解を求めたく，以下のコードを作成しました．
tに値を代入する前はうまく動作しましたが，tをforループで回すとエラーが出てしましました．
どのように改良したらよいでしょうか？

また．16384点の値からRMS値を求めたいのですが、私の書いたコードではおかしい気がします．
どうしたらよいか教えてください．

tををループさせる前のコード（ｔ＝０のときに初期値を与えています）

import sympy as sp
import math
import cv2
import numpy as np
import matplotlib.pyplot as plt
import csv
import os
import glob
import pandas as pd


for i in range(48) :
    Vr=i*0.5+0.5
    print(Vr)
    π = math.pi
    m = 0.0156
    Cn = 6.83879
    D = 0.016
    H = 0.04
    L = 0.3
    ρ = 1.18
    M = (m/(ρ*H*D*L))*(D/H)
    CL = 0.5

    sp.init_printing()
    sp.var('t')
    y = sp.Function('y')(t) #xはtの関数(Function)
    C = (4*π*Cn)/Vr
    F = CL/(2*M)
    eq2 = sp.Eq( sp.diff(y, t ,2)+C*sp.diff(y, t, 1),  F )
    yre = sp.dsolve(eq2, ics={y.subs(t,0):0}) 
    print(yre)

ｔをforを用いて回したときとそのエラー（RMS値のコードも追加）

import sympy as sp
import math
import cv2
import numpy as np
import matplotlib.pyplot as plt
import csv
import os
import glob
import pandas as pd


for i in range(48) :
    Vr = i*0.5 + 0.5
    print(Vr)
    π = math.pi
    m = 0.0156
    Cn = 6.83879
    D = 0.016
    H = 0.04
    L = 0.3
    ρ = 1.18
    M = (m/(ρ*H*D*L))*(D/H)
    CL = 0.5

    for j in range(16384):
        t = j*0.0005
        print(t)
        sp.init_printing()
        y = sp.Function('y')(t) #xはtの関数(Function)
        C = (4*π*Cn)/Vr
        F = CL/(2*M)
        eq2 = sp.Eq( sp.diff(y, t ,2)+C*sp.diff(y, t, 1),  F )
        yre = sp.dsolve(eq2, ics={y.subs(t,0):0}) 
        print(yre)
       
    rms = np.sqrt(np.square(yre).mean(axis=0))
    print(rms)

エラー

ValueError: 
Can't calculate derivative wrt 0.0.

行動規範の内容に同意します

回答1件

ベストアンサー

tへの値を代入後にtで微分しようとしているので、要は定数で微分しようとしてエラーになってるのではないでしょうか。
たとえば下記のように微分後に代入してあげればと思います。（とりあえずすぐにはエラー出ないようですが、長いので最後までは回してません）

ところでこれ、回るは回りますがfor j in range(1, 16384):このfor文は意味をなしていません。
ずっとyreを上書きしているので、最後のj=16383のときのyreだけが残ってrms = np.sqrt(np.square(yre).mean(axis=0))を計算することになるだろうと思いますね。

python3
1import sympy as sp
2import math
3import cv2
4import numpy as np
5import matplotlib.pyplot as plt
6import csv
7import os
8import glob
9import pandas as pd
10
11for i in range(48) :
12    Vr = i*0.5 + 0.5
13    print(Vr)
14    π = math.pi
15    m = 0.0156
16    Cn = 6.83879
17    D = 0.016
18    H = 0.04
19    L = 0.3
20    ρ = 1.18
21    M = (m/(ρ*H*D*L))*(D/H)
22    CL = 0.5
23
24    for j in range(1, 16384):
25        sp.init_printing()
26        sp.var('t')
27        y = sp.Function('y')(t) #xはtの関数(Function)
28        C = (4*π*Cn)/Vr
29        F = CL/(2*M)
30        eq2 = sp.Eq( sp.diff(y, t ,2)+C*sp.diff(y, t, 1),  F )
31        yre = sp.dsolve(eq2, ics={y.subs(t,0):0}) 
32        print(j*0.0005)
33        print(yre)
34        print(yre.subs(t,j*0.0005))
35
36    rms = np.sqrt(np.square(yre).mean(axis=0))
37    print(rms)

投稿2021/07/07 00:51

jeanbiego

総合スコア3966

nonnon615

2021/07/07 01:05

ありがとうございます！＞ところでこれ、回るは回りますがfor j in range(1, 16384):このfor文は意味をなしていません。ずっとyreを上書きしているので、最後のj=16383のときのyreだけが残ってrms = np.sqrt(np.square(yre).mean(axis=0))を計算することになるだろうと思いますね。これを解決し，csvまたはExcelに出力したいのですが，どうしたらよいでしょうか？ずっと上書きされる状態を解決できれば、出力できると思うのですが・・・