【添削のお願い】オイラー法によるシミュレーション

Question

###オイラー法によるシミュレーション
現在、オイラー法によるシミュレーションのコードを書いています。

その結果をエクセルでグラフ化してみました。（下の写真）

すると、周期を細かくしたら、グラフは内に入っていくはずなのですが、
下の画像では、外側にどんどん発散していく形になってしまいます。

どこが、間違っているのでしょうか？

###エクセルで作ったグラフの画像

###問題
ばね定数kのばねの先についている質量mの物体には、
重力mgと、ばねの伸びyに比例したばねの復元力-kyと
速度vに比例した抵抗力-Rv（Rは抵抗係数）が働いている。

【初期条件】
位置　 y0 = 1
初速度　v0 = 0
周期　dt　= 0.0137
m = 1
g = 9.8
k = 22
R = 1.2

###コード

#include  <stdio.h>

int main(void){
   float v1, y1, yh, vh, c1, c2;
   float v0=0;
   float y0=1;
   float k=22;
   float m=1;
   float R=1.2;
   float dt=0.0137;
   float t=0;
   int i=0;
   
   c1 = k/m;
   c2 = R/m;
   
   while(i<300){
       printf("i回目 = %d\n", i);
       printf("t = %f\n", t);
       printf("v0 = %f\n", v0);
       printf("y0 = %f\n", y0);
       printf("\n");
        
        //仮半数先
       yh = y0 + v0 * 0.5 * dt;
       vh = v0 + (-c1*y0 -c2*v0) * 0.5 * dt;
       v1 = v0 + (-c1*yh -c2*vh) * dt;
       y1 = y0 + vh * dt;
       
       v0 = v1;
       y0 = y1;
       i++;
       t = i * dt;
   }
}

以上、よろしくお願いいたします。

Comments

[OHASI]

ご指摘ありがとうございます。

[fana]

> すると、周期を細かくしたら、グラフは内に入っていくはずなのですが、
> 下の画像では、外側にどんどん発散していく形になってしまいます。

「内」とか「外」とか言われてもどういう意味なのかわかりません．
（３つのグラフは単に横軸方向のスケールが異なるだけにしか見えません．
　横軸が　時刻t　ではなく，演算の回数i相当の何か(?)　になってませんか？）

[fana]

少なくとも青とオレンジで描画されたグラフを見るに振幅は減少していっているので

> 速度vに比例した抵抗力-Rv

の効果に依ると思われるエネルギーの減少が起きている（＝「発散」はしていない）ように見えます．

[fana]

あと，

> g = 9.8

コレ使われてないように見えますが…？