2自由度の減衰振動をルンゲクッタ法で計算したいと思い、プログラムを組んでみました。
ただし、実行結果がどうも正しいものになっていないようです。
どなたか詳しい方、教えて下さらないでしょうか?
以下が組んでみたプログラムです。
#include <stdio.h>
#include <math.h>
#define NDIV 500;
static double m1=1.0;
static double m2=1.0;
static double k1=1.0;
static double k2=1.0;
static double c1=0.00;
static double c2=0.00;
// dx/dt = kx
double f(double x1,double v1){
return v1;
}
// dv/dt = kv
double g(double x1,double v1,double x2, double v2){
return -((c1+c2)*v1-c2*v2+(k1+k2)*x1-k2*x2)/m1;
}
double j(double x2, double v2){
return v2;
}
double p(double x1, double v1, double x2, double v2){
return (c2*v1-c2*v2+k2*x1-k2*x2)/m2;
}
int main(void){
double x1; // 変位
double x2;
double v1; // 速度
double v2;
double c1;
double c2; // 減衰
double t; // 時刻
double h; // =dt
double hh; // =dt/2
double kx11,kv11,kx21,kv21;
double kx12,kv12,kx22,kv22;
double kx13,kv13,kx23,kv23;
double kx14,kv14,kx24,kv24;
int i,n=NDIV;
char a[20],*pa;
FILE *fp1,*fp2,*fp3,*fp4,*fp5;
h = 0.1;
hh = h/2;
x1 = 1;
v1 = 0;
x2 = 0;
v2 = 0;
fp1 = fopen("ju-t.txt","w");
fp2 = fopen("ju-x1.txt","w");
fp3 = fopen("ju-v1.txt","w");
fp4 = fopen("ju-x2.txt","w");
fp5 = fopen("ju-v2.txt","w");
if(fp1 == NULL){
printf("can't open/n");
}
else{
printf("open/n");
}
printf(" # t x1 v1 x2 v2\n");
for(i=0;i<=n;i++){
t = h*i;
printf("i,t,x1,v1,x2,v2");
fprintf(fp1, "%.4f\n",t);
fprintf(fp2, "%.4f\n",x1);
fprintf(fp3, "%.4f\n",v1);
fprintf(fp4, "%.4f\n",x2);
fprintf(fp5, "%.4f\n",v2);
kx11 = f(x1,v1);
kv11 = g(x1,v1,x2,v2);
kx21 = j(x2,v2);
kv21 = p(x1,v1,x2,v2);
kx12 = f(x1+hh*kx11,v1+hh*kv11);
kv12 = g(x1+hh*kx11,v1+hh*kv11,x2+hh*kx21,v2+hh*kv21);
kx22 = j(x2+hh*kx21,v2+hh*kv21);
kv22 = p(x1+hh*kx11,v1+hh*kv11,x2+hh*kx21,v2+hh*kv21);
kx13 = f(x1+hh*kx12,v1+hh*kv12);
kv13 = g(x1+hh*kx12,v1+hh*kv12,x2+hh*kx22,v2+hh*kv22);
kx23 = j(x2+hh*kx22,v2+hh*kv22);
kv23 = p(x1+hh*kx12,v1+hh*kv12,x2+hh*kx22,v2+hh*kv22);
kx14 = f(x1+h*kx13,v1+h*kv13);
kv14 = g(x1+h*kx13,v1+h*kv13,x2+h*kx23,v2+h*kv23);
kx24 = j(x2+h*kx23,v2+h*kv23);
kv24 = p(x1+h*kx13,v1+h*kv13,x2+h*kx23,v2+h*kv23);
x1 = x1 + h*(kx11+2*kx12+2*kx13+kx14)/6;
v1 = v1 + h*(kv11+2*kv12+2*kv13+kv14)/6;
x2 = x2 + h*(kx21+2*kx22+2*kx23+kx24)/6;
v2 = v2 + h*(kv21+2*kv22+2*kv23+kv24)/6;
}
fclose(fp1);
fclose(fp2);
fclose(fp3);
fclose(fp4);
fclose(fp5);
return 0;
}
A 回答 (2件)
- 最新から表示
- 回答順に表示
No.2
- 回答日時:
「過去類似の研究を行っていた先輩が残してくれた実行結果と異なっていた」としても, パラメータが違っていたら当然違う形になりえるわけだからそれが直接「実行結果が間違っている」ことにはつながらないのでは?
で
http://oshiete.goo.ne.jp/qa/8642263.html
は放置ですか?
No.1
- 回答日時:
どういう微分方程式を解こうとしているのですか?
「実行結果がどうも正しいものになっていないようです」の根拠は?
この回答への補足
解こうとしている微分方程式は
m1*x1" = -(c1+c2)x1'+c2*x2'-(k1+k2)*x1+k2*x2
m2*x2" = c2*x1'-c2*x2'+k2*x1-k2*x2
の連立微分方程式です。
'は一階時間微分、"は二階時間微分のつもりです。
また、kはバネ係数、cは減衰係数です。
自分は大学生の研究室所属なのですが、過去類似の研究を行っていた先輩が
残してくれた実行結果と異なっていたので、正しくないのではと判断しました。
お探しのQ&Aが見つからない時は、教えて!gooで質問しましょう!
似たような質問が見つかりました
- C言語・C++・C# ある線が円の範囲に入っているかの計算 1 2022/12/07 16:14
- C言語・C++・C# プログラミングの授業の課題です 1 2023/01/17 22:15
- C言語・C++・C# C言語の課題が出たのですが自力でやっても分かりませんでした。 要素数がnであるint型の配列v2の並 3 2022/11/19 17:41
- C言語・C++・C# バイナリファイルをコピーするのにかかる時間を測りたいのですが実行するとFatel error:gli 2 2022/11/03 01:10
- Visual Basic(VBA) VBAプログラミング 2 2022/11/27 12:07
- Visual Basic(VBA) VBAプログラミング 2 2022/11/27 12:13
- 物理学 写真の図は単振動の動きを段階的に表したものです。 (加速度=a、力=kx、ばね定数=k、物体の質量= 8 2022/08/24 23:39
- 物理学 ばねの単振動による復元力はf=-kxと表され、フックの大きさによって表される弾性力はf=kxと表され 4 2023/05/27 00:13
- C言語・C++・C# c言語でユーザ関数を利用して複素数のべき乗と絶対値の数列を計算するプログラムが作りたいです。 3 2023/01/29 22:13
- その他(プログラミング・Web制作) Pythonにおける物理のシミュレーションでの単位変換について 2 2023/06/02 17:11
このQ&Aを見た人はこんなQ&Aも見ています
関連するカテゴリからQ&Aを探す
おすすめ情報
このQ&Aを見た人がよく見るQ&A
デイリーランキングこのカテゴリの人気デイリーQ&Aランキング
-
プログラムでの数字につく”f”の...
-
doubleの変数にintとintの割り...
-
c言語で、繰り返し文の中で、0....
-
C言語を実行すると-infが出てき...
-
float型とdouble型の変数の違い...
-
C言語の型による処理速度の違い
-
C言語 関数プロトタイプ宣言の...
-
C言語で-23乗を取り扱うには
-
方程式を2分法を用いて解くプロ...
-
EXE1→DLL→EXE2数値を受け渡す方法
-
C言語でdouble型の小数点の引き...
-
浮動小数点の定数
-
difftime()について
-
(C,C++言語)関数の引数は自動キ...
-
mallocで動的確保後、値が変わる
-
DWORDの警告
-
数値を指数部と仮数部に分離したい
-
c言語の問題
-
C言語でプログラミングしまし...
-
C++で外積
マンスリーランキングこのカテゴリの人気マンスリーQ&Aランキング
-
プログラムでの数字につく”f”の...
-
doubleの変数にintとintの割り...
-
float型とdouble型の変数の違い...
-
c言語で、繰り返し文の中で、0....
-
C言語を実行すると-infが出てき...
-
C言語 関数プロトタイプ宣言の...
-
C 開放してるのにエラー(doubl...
-
C言語の型による処理速度の違い
-
至急です! マクロ定義で #defi...
-
関数におけるif文とreturn文に...
-
c言語のプログラミングについて...
-
2分法で方程式の複数の解を自...
-
-1.#IND00と出てしまうのですが...
-
doubleは常に%lfとするべきなのか
-
C言語のpow関数の不具合
-
C言語で-23乗を取り扱うには
-
C言語で台形公式を使った二重積...
-
Cで3乗根を求める方法
-
sin(x)の近似について
-
2次方程式の解を求めるプログ...
おすすめ情報