※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 の不定積分がわ...
-
∫1/√x dx 積分せよ 教えて下さい
-
項の右端につく縦棒の意味を教...
-
フーリエ級数の問題で、f(x)は...
-
【数学Ⅱ・Ⅲ】微分の問題
-
積分 Xの-2乗を積分するとどう...
-
x−1分の2の微分の仕方を教えて...
-
∫e^cos(x) dx の計算
-
1/X^2の積分ってlogX^2ですか?
-
理経数学プラチカ55番(2)で・(d...
-
ある積分の問題。∫1/√(x^2+A) =...
-
定積分∫[3→0]|x^2-4|dxの答え...
-
微積分 dの意味
-
sinθの軌跡距離は?
-
1/(x*(x-1)^(1/2))の積分について
-
微分方程式
-
虚数「i」の無限大への極限
-
これが解けません。。。
マンスリーランキングこのカテゴリの人気マンスリーQ&Aランキング
-
積分で1/x^2 はどうなるのでし...
-
e^-2xの積分
-
フーリエ級数の問題で、f(x)は...
-
∫1/(x^2+1)^2 の不定積分がわ...
-
積分 Xの-2乗を積分するとどう...
-
∫1/√x dx 積分せよ 教えて下さい
-
項の右端につく縦棒の意味を教...
-
1/X^2の積分ってlogX^2ですか?
-
何が違いますか?
-
微積分 dの意味
-
∫e^cos(x) dx の計算
-
【数学Ⅱ・Ⅲ】微分の問題
-
2次微分の変数変換
-
x/(a^2+x^2)の積分について
-
exp(-ax^2)*cosx の証明
-
∫r/(a^2+r^2)^3/2drの計算の解...
-
x^2 * exp(x^2) dxの不定積分
-
x−1分の2の微分の仕方を教えて...
-
フーリエ変換の問題について
-
台形の任意の高さにおける上辺...
おすすめ情報