![](http://oshiete.xgoo.jp/images/v2/pc/qa/question_title.png?e8efa67)
※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
![「単振り子とルンゲ・タック法」の質問画像](http://oshiete.xgoo.jp/_/bucket/oshietegoo/images/media/f/542948080_62d02ee8972dc/M.png)
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ランキング
-
e^-2xの積分
-
積分で1/x^2 はどうなるのでし...
-
∫1/(x^2+1)^2 の不定積分がわ...
-
∫1/√x dx 積分せよ 教えて下さい
-
フーリエ級数の問題で、f(x)は...
-
1/X^2の積分ってlogX^2ですか?
-
ガウス積分の計算がうまくいき...
-
積分 Xの-2乗を積分するとどう...
-
2次微分の変数変換
-
台形の任意の高さにおける上辺...
-
【数学Ⅱ・Ⅲ】微分の問題
-
xe^xsinx これを部分積分するに...
-
∮a^xdxこれを公式的に導いてほ...
-
微積分 dの意味
-
x/(a^2+x^2)の積分について
-
項の右端につく縦棒の意味を教...
-
exp(-ax^2)*cosx の証明
-
(dy/dx)+y=xの微分方程式はどの...
-
数学の微分の問題
-
定積分の問題です。
マンスリーランキングこのカテゴリの人気マンスリー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)は、階乗からど...
-
フーリエ変換の問題について
-
x/(a^2+x^2)の積分について
-
(dy/dx)+y=xの微分方程式はどの...
-
x−1分の2の微分の仕方を教えて...
-
虚数「i」の無限大への極限
-
微分方程式の解き方
おすすめ情報