C言語 ルンゲクッタ法２階微分の計算

Question

ご覧いただきありがとうございます。

ルンゲクッタ法を用いて

dq/dt=p
dp/dt=-q
p0=1
q0=0 

を解きたいです。
現在、ウェブサイトでルンゲクッタ法の解き方を自分で調べて書いてみたのですが、明らかにおかしな値（sin,cos の値になっていない）になっています。
参考にしたところに書いてあったのは

d^2y/dx^2=f(x,y,y')
を解くのに
y'=p
とおいて
dp1=h*f(x,y,p)
dy1=h*p
dp2=h*f(x+h/2,y+dy1/2,p+dp1/2)
dy2=h*(p+dp1/2)
dp3=h*f(x+h/2,y+dy2/2,p+dp2/2)
dy3=h*f(p+dp2/2)
dp4=h*f(x+h,y+dy3,p+dp3)
dy4=h*(p+dp3)

で私の場合
y=q,x=t,y'=p,f(x,y,y')=-q
と置き換えました。
こう置き換えた場合、p,qがtの関数であるはずなのにうまくtでパラメータ表示できていないのですが、どう改善すればいいのでしょうか。

double runge(double q0, double p0, double t0, double tn, int n)
{
    int i;
    double q, p, t, h, dq1, dq2, dq3, dq4, dp1, dp2, dp3, dp4;
    q = q0;
    p = p0;
    t = t0;
    h = (tn - t0) /n;
    
    
    for ( i=1; i <= n ; i++){
      
        dp1 = -q;
        dq1 = p;
        dp2 = -(q+dq1/2);
        dq2 = p + dp1/2;
        dp3 = -(q + dq2/2);
        dq3 = p+dp2/2;
        dp4 = -(q + dq3);
        dq4 = p+ dp3;
        q += (dq1 + 2 * dq2 + 2 * dq3 + dq4)*(h/6.0); 
　　　t = t0 + i*h;
        printf("q(%f)=%f\n", t, q);
 　　 p += (dp1 + 2 * dp2 + 2 * dp3 + dp4)*(h/6.0); 
        printf("p(%f)=%f\n", t, p);

Comments

[RiaFeed]

forの外でtに値を代入しているのにforのすぐ下でまたtに違う値を代入していますね。
多分ループ1回分ずれてるんじゃないかな？

[hoshi-takanori]

Objective-C は C とは別の言語なので、C のタグを使ってください。 https://teratail.com/tags/C

[fana]

> 明らかにおかしな値になっています

って言うけど，じゃあ「正しい値」っていくつなんですか？　っていうのを示さないと伝わらないのでは．

で，それはそれとして，pを更新する処理が無いように見えますが．

[yusuke.]

ご指摘ありがとうございます。
pについてはコピーミスです。

[fana]

区間幅 h が各勾配値の計算にも現れねばならないと思うのですが，どうなのでしょう？

Answer 1

正直，自信はありませんが，以下の考えで実装してみたところ，それっぽい値が得られた感じです．

---

[考え]

（いろいろ省略せずに）数式を以下のように書いてみましょう．

P'(t) = -Q( t, P(t) )
Q'(t) = P( t, Q(t) )

すなわち，右辺の関数の引数に P(t) や Q(t) が含まれている形に．

で，ルンゲクッタ法について考えてみます．
関数 P(t) の側を例にすると，次の時刻での値 P(t+dt) を求めるために，

1. 時刻 t における勾配 P'(t) （: これは数式で与えられている）
2. 1.の勾配を用いて推定した区間中点での勾配 P'(t+dt/2)
3. 2.の勾配を用いて推定した区間中点での勾配 P'(t+dt/2)
4. 3.の勾配を用いて推定した区間末尾での勾配の推定値 P'( t+dt )

を用いるわけで，これらの計算をどうするのか？　っていうのが問題点でしょう．

例えば，2.について考えてみる（3.や4.も話は同様なので）と，求める P'(t+dt/2) というのは，最初に書いたように，
P'(t+dt/2) = -Q( t+dt/2, P(t+dt/2) ) ですから，

• ここの右辺にある P(t+dt/2) の値はどうすればよいのか？
• そしてそれをどう使うのか？

を考える必要がありそうですね．
ルンゲクッタ法のアルゴリズムの意味に沿って考えれば，このステップ(2.)というのは，
【1.の勾配を用いて，区間中点での勾配をオイラー法的に求める】ですから，その値については
P(t) + 0.5*dt + (1.で求めた勾配)
とすれば良さそうです．

では，この値の用途はどうするのか？
すなわち，この値をどうにかして用いて Q( t+dt/2, この値 ) の推定値を作る必要があるのだけど，どうするのか？っていう．
当然(?)，ここでの Q(t+dt/2) もオイラー法的に｛ Q(t), P(t) と，この P'(t+dt/2) 推定値 ｝から作ればよいのだと思いますが，しかしながら素のオイラー法を考えると P'(t+dt/2) 推定値の出番がなくなってしまう^^
→なので，ここを「修正オイラー法な感じにすることで，こいつを使ってやるぜ！」とすればよいのではないかな，と考えました．すなわち，
P'(t+dt/2) = -Q( t+dt/2, P(t+dt/2) ) の右辺にある Q( t+dt/2, P(t+dt/2) ) の値を
= Q(t) + 0.5*dt * ( Q'(t) + Q'(t+dt)推定値 )*0.5 として２つの勾配値から求めることとし，ここで Q'(t+dt)推定値 = P(t+dt/2)推定値 として使用しています．

---

[コード]

int main()
{
	double t = 0;	//※tはprintf()にしか使ってないです．
	double P_t = 1;	//初期値 P(t=0) = 1
	double Q_t = 0;	//初期値 Q(t=0) = 0

	//適当な刻み幅でやってみた
	const double PI = acos(-1.0);
	const double dt = PI/180.0;

	for( int i=0; i<360; ++i )
	{
		//P(t+dt)をルンゲクッタ法で求む
		double P_next;
		{//P_next
			//tでの勾配
			double dP_t = -Q_t;	//P'(t) = -Q(t)

			//中点 t+dt/2 での勾配の推定値２つ
			double dP_mid[2];	// P'(中点) = -Q( 中点, P(中点)推定値 )
			{//dP_mid[0]
				//dP_t を用いて推定した P(中点)
				double P_mid_est = P_t + 0.5*dt*dP_t;

				//Q(t)に関する２つの勾配値
				//  Q'(t) = P(t)
				//  Q'(中点) = P(中点)
				//の２つから Q(t+dt/2) を推定（修正オイラー法な感じで）
				double Q_mid_est = Q_t  +  0.5 * 0.5*dt * ( P_t + P_mid_est );
				//
				dP_mid[0] = -Q_mid_est;
			}
			{//dP_mid[1]
				//dP_mid[0] を用いて推定した P(中点)
				double P_mid_est = P_t + 0.5*dt*dP_mid[0];
				//(dP_mid[0] と同様)
				double Q_mid_est = Q_t  +  0.5 * 0.5*dt * ( P_t + P_mid_est );
				dP_mid[1] = -Q_mid_est;
			}

			//区間末尾 t+dt での勾配推定値
			double dP_next;	//P'(末尾) = -Q( 末尾, P(末尾)の推定値 )
			{
				//dP_mid[1] を用いて推定した P(末尾)
				double P_next_est = P_t + dt*dP_mid[1];
				//
				double Q_next_est = Q_t  +  0.5 * dt * ( P_t + P_next_est );
				dP_next = -Q_next_est;
			}

			//
			P_next = P_t + dt*( dP_t + 2*( dP_mid[0] + dP_mid[1] ) + dP_next )/6.0;
		}

		//Q(t+dt)側も同様
		double Q_next;
		{//Q_next
			//tでの勾配
			double dQ_t = P_t;	//Q'(t) = P(t)

			//中点 t+dt/2 での勾配の推定値２つ
			double dQ_mid[2];	// Q'(中点) = P( 中点, Q(中点)推定値 )
			{//dQ_mid[0]
				double Q_mid_est = Q_t + 0.5*dt*dQ_t;

				//P(t)に関する２つの勾配値
				//  P'(t) = -Q(t)
				//  P'(中点) = -Q(中点)
				//の２つから P(t+dt/2) を推定（修正オイラー法な感じで）
				double P_mid_est = P_t  +  0.5 * 0.5*dt * -( Q_t + Q_mid_est );
				dQ_mid[0] = P_mid_est;
			}
			{//dQ_mid[1]
				double Q_mid_est = Q_t + 0.5*dt*dQ_mid[0];
				double P_mid_est = P_t  +  0.5 * 0.5*dt * -( Q_t + Q_mid_est );
				dQ_mid[1] = P_mid_est;
			}

			//区間末尾 t+dt での勾配推定値
			double dQ_next;	//Q'(末尾) = P( 末尾, Q(末尾)の推定値 )
			{
				double Q_next_est = Q_t + dt*dQ_mid[1];
				double P_next_est = P_t  +  0.5 * dt * -( Q_t + Q_next_est );
				dQ_next = P_next_est;
			}

			//
			Q_next = Q_t + dt*( dQ_t + 2*( dQ_mid[0] + dQ_mid[1] ) + dQ_next )/6.0;
		}

		//更新
		P_t = P_next;
		Q_t = Q_next;
		t += dt;
		printf( "%f, %f, %f\n", t, P_t, Q_t );
	}

	return 0;
}