※VBAカテにするかどうか迷ったのですが、微分方程式の数値解の話ですのでこちらで質問させてください。
Excel VBA で単振り子の厳密解をルンゲ・タック法で解いてみたのですが、初期条件とした与えた角度(単振動で言えば振幅)が時間が経過すると共に大きくなっていきます。どこに問題があるのでしょうか?
ここの画像は小さいので(いい加減改善してほしい)
http://imepic.jp/20220714/862800
でご確認願います。
Sub SiPendulum()
Worksheets("単振り子").Activate
Dim K1 As Double, K2 As Double, K3 As Double, K4 As Double
Dim L1 As Double, L2 As Double, L3 As Double, L4 As Double
Dim t As Double, x As Double, v As Double
Dim h As Double, Dx As Double, Dv As Double
Dim i As Long
Cells(5, 5).Value = "t"
Cells(5, 6).Value = "θ"
Cells(5, 7).Value = "ω"
' d^2θ/dt^2 = -(g/L)sin(θ) ------> dω/dt = -(g/L)sin(θ).
' dθ/dt = ω
' 変数名を変える θ→x ω→v
' d^2x/dt^2 = -(g/L)sin(x) ------> dv/dt = -(g/L)sin(x).
' dx/dt = v
h = 0.05
t = 0
x = Cells(5, 11).Value '初期位相
v = 0
Cells(6, 5).Value = t
Cells(6, 6).Value = x
Cells(6, 7).Value = v
For i = 1 To 81
K1 = h * Px_t(t, x, v) 'dx/dt
L1 = h * Pv_t(t, x, v) 'dv/dt
K2 = h * Px_t(t + h / 2, x + K1 / 2, v + L1 / 2)
L2 = h * Pv_t(t + h / 2, x + L1 / 2, v + L1 / 2)
K3 = h * Px_t(t + h / 2, x + K2 / 2, v + L1 / 2)
L3 = h * Pv_t(t + h / 2, x + L2 / 2, v + L1 / 2)
K4 = h * Px_t(t + h, x + K3, v + L3)
L4 = h * Pv_t(t + h, x + L3, v + L3)
Dx = (K1 + 2 * K2 + 2 * K3 + K4) / 6
Dv = (L1 + 2 * L2 + 2 * L3 + L4) / 6
t = t + h
x = x + Dx
v = v + Dv
Cells(6 + i, 5).Value = t
Cells(6 + i, 6).Value = x
Cells(6 + i, 7).Value = v
Next i
End Sub
'dx/dt = v
Function Px_t(t, x, v)
Px_t = v
End Function
'dv/dt = -(g/L)sin(x)
Function Pv_t(t, x, v)
Worksheets("単振り子").Activate
Dim GL As Double
GL = Cells(6, 4).Value
Pv_t = -GL * Sin(x)
End Function
No.1ベストアンサー
- 回答日時:
古典的な 4段 4次のルンゲ・クッタ法だと思うけど, L2, K3, L3, L4 の計算は間違ってない? 例えば, L2 の式
L2 = h * Pv_t(t + h / 2, x + L1 / 2, v + L1 / 2)
は x+L1/2 でいいの?
あと t の関数 x(t) について dx/dt が陽に現れない微分方程式
d^2x/dt^2 = f(t, x)
だと同じ 4段 4次でももうちょっと精度を上げることができたはず.
> L2 = h * Pv_t(t + h / 2, x + L1 / 2, v + L1 / 2)
> は x+L1/2 でいいの?
あちゃ~、お恥ずかしい(笑)。
> あと t の関数 x(t) について dx/dt が陽に現れない微分方程式
> d^2x/dt^2 = f(t, x)
> だと同じ 4段 4次でももうちょっと精度を上げることができたはず.
知りませんでした。感謝いたします。
お探しのQ&Aが見つからない時は、教えて!gooで質問しましょう!
似たような質問が見つかりました
- Visual Basic(VBA) VBAプログラミング 2 2022/11/27 12:07
- Visual Basic(VBA) VBAプログラミング 2 2022/11/27 12:13
- Visual Basic(VBA) VBA 別ブックからの転記の高速化について VBA 別ブックからの転記の高速化についてご教授下さい。 19 2022/07/26 13:07
- Visual Basic(VBA) いつもお世話になっております、VBAで教えて頂きたいのですが 2 2022/05/05 22:20
- Excel(エクセル) VBAの指示の内容 昨日こちらでご教示頂いたのですが初心者な為、一つ一つの指示が何をやっているのかわ 2 2022/10/25 18:08
- Visual Basic(VBA) VBA処理追加 こちらでご教示頂いたのですが回答完了させてしまいましたのでこちらからまた質問させてく 2 2022/10/27 09:57
- Visual Basic(VBA) VBAで質問ですが、皆さんはどの様に導き出しているのでしょうか? 6 2022/05/03 21:53
- Visual Basic(VBA) vbaの計算 if elseと範囲について 6 2022/11/26 01:49
- Visual Basic(VBA) 別シートから年齢別の件数をカウントしたいの続き 5 2023/01/24 00:16
- Visual Basic(VBA) vbaを早くしたい 5 2022/09/09 10:58
デイリーランキングこのカテゴリの人気デイリーQ&Aランキング
-
積分で1/x^2 はどうなるのでし...
-
e^-2xの積分
-
∫1/√x dx 積分せよ 教えて下さい
-
∫1/(x^2+1)^2 の不定積分がわ...
-
積分 Xの-2乗を積分するとどう...
-
部分積分の話ですが赤く印をつ...
-
∫e^cos(x) dx の計算
-
2次微分の変数変換
-
フーリエ級数の問題で、f(x)は...
-
項の右端につく縦棒の意味を教...
-
単振り子とルンゲ・タック法
-
(1) dy/dx=1-y^2 (2)dy/dx=y-y^...
-
微分の記号
-
ガウス関数のフーリエ変換 証明
-
【数学Ⅱ・Ⅲ】微分の問題
-
偶関数と奇関数の見分け方
-
媒介変数(t)の2回微分について...
-
導関数の記号 dy/dx の意味は?
-
d^2x/dt^2=-ω^2x-2γ(dx/dt) (ω=...
-
確率密度関数をf(x)=1-|x-1|と...
マンスリーランキングこのカテゴリの人気マンスリーQ&Aランキング
-
積分で1/x^2 はどうなるのでし...
-
e^-2xの積分
-
∫1/(x^2+1)^2 の不定積分がわ...
-
∫1/√x dx 積分せよ 教えて下さい
-
フーリエ級数の問題で、f(x)は...
-
積分 Xの-2乗を積分するとどう...
-
項の右端につく縦棒の意味を教...
-
1/X^2の積分ってlogX^2ですか?
-
∫e^cos(x) dx の計算
-
微積分 dの意味
-
【数学Ⅱ・Ⅲ】微分の問題
-
これはわかる
-
x−1分の2の微分の仕方を教えて...
-
フーリエ変換の問題について
-
1階微分方程式の解析解に関する...
-
x/(a^2+x^2)の積分について
-
2次微分の変数変換
-
(dy/dx)+y=xの微分方程式はどの...
-
x^2 * exp(x^2) dxの不定積分
-
e^-1/Tの積分
おすすめ情報