コンピュータ(プログラミング)のカテゴリに投稿しようかとも考えましたが、物理学に関する数値計算ということなので、物理学のカテゴリに投稿しました。
単振動や減衰振動、そして強制振動などを数値解析したいと思い、ルンゲクッタ法を使いシミュレートしてみようとしています。ルンゲクッタ法という方法を全く知らなかったので、インターネットや図書館で調べたのですが、どうしても分からないことろがあるので質問することにしました。
書籍やネットの情報を参考にしながら、単振動の場合を数値解析してみました(C言語を使って)。この単振動はうまくできたのですが、どうしても、その先の、減衰振動の数値解析がうまくいかないので、困っています。
----
#include <stdio.h>
double f1(double t,double x,double v);
double f2(double t,double x,double v);
int main()
{
double t,x,v,dt,tmax;
double k1[2],k2[2],k3[2],k4[2];
x=1.0; //位置の初期値
v=0.0; //速度の初期値
dt=0.01; //刻み幅
tmax=500.0; //繰り返し最大回数
FILE *output;
output=fopen("output.dat","w");
for(t=0;t<tmax;t+=dt) {
k1[0]=dt*f1(t,x,v);
k1[1]=dt*f2(t,x,v);
k2[0]=dt*f1(t+dt/2.0,x+k1[0]/2.0,v+k1[1]/2.0);
k2[1]=dt*f2(t+dt/2.0,x+k1[0]/2.0,v+k1[1]/2.0);
k3[0]=dt*f1(t+dt/2.0,x+k2[0]/2.0,v+k2[1]/2.0);
k3[1]=dt*f2(t+dt/2.0,x+k2[0]/2.0,v+k2[1]/2.0);
k4[0]=dt*f1(t+dt,x+k3[0],v+k3[1]);
k4[1]=dt*f2(t+dt,x+k3[0],v+k3[1]);
x=x+(k1[0]+2.0*k2[0]+2.0*k3[0]+k4[0])/6.0;
v=v+(k1[1]+2.0*k2[1]+2.0*k3[1]+k4[1])/6.0;
fprintf(output,"%f %f %f\n",t,x,v);
}
fclose(output);
return 0;
}
double f1(double t,double x,double v)
{
return v;
}
double f2(double t,double x,double v)
{
return (-x);
}
----
このソースは単振動のもので、減衰振動のときは、最後の方の
----
double f1(double t,double x,double v)
{
return v;
}
double f2(double t,double x,double v)
{
return (-x);
----
の部分が変わるのだと思うのですが、よく分かりません。
減衰振動は、mx"=-kx-rx'で表されるので、この式を変形(?)したものが、入るのかな、という予想程度にしか分かりません。
ソースをどのように変えれば減衰振動を解析できるのでしょうか。
どなたか詳しい方、教えてください。お願いします。
A 回答 (2件)
- 最新から表示
- 回答順に表示
No.2
- 回答日時:
そんだけわかっているんだから, そのまま書けばいいだけだろうに....
単振動の場合は mx'' = -kx で, v = dx/dt を導入することで
dx/dt = v, m dv/dt = -kx
という連立 1階微分方程式にしたんですよね? 減衰振動の mx'' = -kx - rx' も同じように連立 1階微分方程式にするだけです. つまり, こちらでも v = dx/dt を導入して
dx/dt = f1(t, x, v), dv/dt = f2(t, x, v)
の形にするだけ.
お探しのQ&Aが見つからない時は、教えて!gooで質問しましょう!
似たような質問が見つかりました
- 数学 単振り子とルンゲ・タック法 1 2022/07/15 00:05
- C言語・C++・C# Cのdoubleの浮動小数点表示について 3 2023/04/17 13:14
- C言語・C++・C# C 言語の Gauss Jordan 法について 2 2022/12/28 11:16
- C言語・C++・C# c言語でユーザ関数を利用して複素数のべき乗と絶対値の数列を計算するプログラムが作りたいです。 3 2023/01/29 22:13
- C言語・C++・C# C言語のマクローリン展開ローラン展開のコードについて 3 2022/12/15 14:45
- C言語・C++・C# C++のcinの動作 5 2023/02/26 00:13
- C言語・C++・C# ある線が円の範囲に入っているかの計算 1 2022/12/07 16:14
- その他(プログラミング・Web制作) 物理の斜方投射のシミュレーションにおける位置や速度の単位について 4 2023/05/31 09:50
- C言語・C++・C# C言語初心者 構造体 課題について 1 2023/03/10 19:30
- C言語・C++・C# バイナリファイルをコピーするのにかかる時間を測りたいのですが実行するとFatel error:gli 2 2022/11/03 01:10
このQ&Aを見た人はこんなQ&Aも見ています
おすすめ情報
このQ&Aを見た人がよく見るQ&A
デイリーランキングこのカテゴリの人気デイリーQ&Aランキング
-
d^2r/dt^2の意味
-
電流の時間微分、電圧の時間微分
-
雨滴の運動質量が変化する落体...
-
力学の雨滴の落下の問題です
-
EXCEL上の数字を自動で振り分け...
-
質量流量の記号「・ の読み方を...
-
力学について質問です。 1.棒の...
-
ノイマン境界条件の球の拡散方...
-
運動方程式の動径方向と方位角方向
-
ルンゲクッタ法で数値解析(C...
-
運動方程式の成分表示
-
ポテンシャルエネルギーから力...
-
伝達関数を求めることができる...
-
∫1/x√(x+1)dx を教えてください!
-
蒸発速度
-
機械力学の問題です!!!
-
物理学の問題です。 図に示すよ...
-
物理の計算で m×dv/dt×v=d/dt{...
-
物理学 おそらく剛体の問題です
-
v^2-v0^2=2ax 今日この式を習っ...
マンスリーランキングこのカテゴリの人気マンスリーQ&Aランキング
-
電流の時間微分、電圧の時間微分
-
d^2r/dt^2の意味
-
ラプラス方程式
-
Debug.Printで表示される内容を...
-
極座標の運動方程式の計算の間...
-
質量流量の記号「・ の読み方を...
-
微分積分のdの意味
-
最後のdv/dtは何でしょうか。
-
加速度 a=dv/dt = (d^2 x) /dt^2
-
雨滴の運動質量が変化する落体...
-
力学について質問です。 1.棒の...
-
運動方程式を求めてください
-
機械力学の問題です!!!
-
微分記号“d”について
-
v^2-v0^2=2ax 今日この式を習っ...
-
運動方程式の微分積分の計算
-
運動方程式の成分表示
-
ブックオフは潰れないよね? キ...
-
EXCEL上の数字を自動で振り分け...
-
難しくてわからないので教えて...
おすすめ情報