※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^2+1)^2 の不定積分がわ...
-
フーリエ級数の問題で、f(x)は...
-
∫1/√x dx 積分せよ 教えて下さい
-
微分方程式 について d²y/dx² ...
-
項の右端につく縦棒の意味を教...
-
積分 Xの-2乗を積分するとどう...
-
台形の任意の高さにおける上辺...
-
1/X^2の積分ってlogX^2ですか?
-
【数学Ⅱ・Ⅲ】微分の問題
-
フーリエ変換の問題について
-
虚数「i」の無限大への極限
-
微積分 dの意味
-
∫e^cos(x) dx の計算
-
大学数学「空間ベクトル」の問...
-
広義積分の値
-
x/(a^2+x^2)の積分について
-
∫r/(a^2+r^2)^3/2drの計算の解...
-
定積分∫[3→0]|x^2-4|dxの答え...
マンスリーランキングこのカテゴリの人気マンスリーQ&Aランキング
-
積分で1/x^2 はどうなるのでし...
-
e^-2xの積分
-
∫1/(x^2+1)^2 の不定積分がわ...
-
∫1/√x dx 積分せよ 教えて下さい
-
フーリエ級数の問題で、f(x)は...
-
積分 Xの-2乗を積分するとどう...
-
項の右端につく縦棒の意味を教...
-
∫e^cos(x) dx の計算
-
1/X^2の積分ってlogX^2ですか?
-
微積分 dの意味
-
exp(-ax^2)*cosx の証明
-
2次微分の変数変換
-
【数学Ⅱ・Ⅲ】微分の問題
-
フーリエ変換の問題について
-
ガンマ関数Γ(x)は、階乗からど...
-
(dy/dx)+y=xの微分方程式はどの...
-
x/(a^2+x^2)の積分について
-
x−1分の2の微分の仕方を教えて...
-
虚数「i」の無限大への極限
-
微分方程式の解き方
おすすめ情報