• 締切済み

ルンゲクッタ法による岩石の軌跡を求める

ルンゲクッタ法により、地球の周りの岩石がどのような輪の軌跡で円運動しているのかを求める勉強をしています。 が、プログラムは作成したものの、結果がまったくありえない数字になってしまったのです… バックグラウンドとしては 岩石は地球から500kmの付近を10km/sで移動している。 微分方程式はf(x,v)=v(t) g(x,v)=-GMx(t)/r*r*r これをy成分でも行なう。 そしてxとyの座標をプロットする。 Gは重力定数、Mは地球の重力を表しています。 このxとyの座標が時間(t)によってどのような軌跡を描いて移動するのかがわかればいいのですが…予測される綺麗な円運動とは程遠い結果になってしまいました… どなたか、「この分野だったら任せろ!」「このHPなら詳しくのってるよ~」という方はいらっしゃいませんでしょうか? 私自身が検索しても、振り子の単振動などしか見つけることができませんでした。似たような例が載っているHPなどありましたら是非教えてください。お願いいたします。 長文失礼いたしました。

  • zamax
  • お礼率11% (1/9)

みんなの回答

回答No.3

 この設定での解析解を求めることは、気の利いた微分方程式の教科書には載っていると思いますが、数値解である必要はあるのでしょうか?  また、逆二乗則の運動方程式を解くのならルンゲクッタより適切な手法があると思いますが、ルンゲクッタでなければならないのでしょうか?

回答No.2

#1です。 申し訳ありません。 見慣れない表記だったもので、勘違いをしました。 微分方程式は、間違いではなく、私の読み違いのようですね。 大変失礼したしました。お詫び申し上げます。プログラムをテストして、適切な刻み幅を選択してみてください。

回答No.1

Runge-Kutta法の問題ではなく、物理の問題になるかと思います。 ・まず、岩石の軌跡はきれいな円運動にはなりません。一般に楕円になります(初期値によっては、きれいな円運動にもなり得ますが)。 ・次に、微分方程式が間違っていますね。万有引力の法則をもう一度見直して下さい。 それでもおかしい場合は、プログラムそのものを疑いましょう。簡単な微分方程式の例題(ご自分で解けるもので十分です)で何回かテストして見て下さい。まあ、本来ならテストは終わっている筈ですけれど。 Runge-Kutta法の話も簡単にしておきましょう。 重力可の運動は、微分方程式に特に困難な点がないため、刻み幅(h)の選択も神経を使う必要はなく、適当に小さい値を固定したまま使って問題ないでしょう。適当な刻み幅で始めてみて、解が安定するまで刻み幅を半分に減らし続けていけば、まあまあ悪くない刻み幅が得られます。今勉強されている範囲であれば、それで十分でしょう。 ご不明の点は、補足でお願いします。

zamax
質問者

お礼

ありがとうございました! 値を全体的に小さくしたらちゃんと円(楕円)になりましたよ! 次に離心率を求めるんですが、これがまた苦労してます。

関連するQ&A

  • 軌跡のy座標最大

    ※→rはrベクトルの事。[ ]内は添字。 →r=(v[x0]t+x[0],g/2t^2+v[y0]+y[0])で与えられる時、速度と加速度を求めよ。 但し、v[x0],x[0],g,v[y0],y[0]は定数。 (1) v[x]=v[x0] v[y]=-gt+v[y0] a[x]=0 a[y]=-g (2) 位置x,yの軌跡の式y=f(x)を求めよ。 y=-g/2(x-x[0]/v[x0])^2+v[y0](x-[x0]/v[x0])+y[0] (3) 軌跡のy座標が最大となるxを求めよ、という問題が分かりません。 お願いします。

  • 軌跡と方程式

    『aの値が変化するとき、y=x~2-4ax+1の頂点Pの軌跡を求めよ。』 という問題なのですが、aの座標を(u,v)とし、y=x~2-4ax+1を平方完成した結果のx座標とy座標に代入してみたのですが、どうも違うようです。 解答によるとaの座標は(2a,-4a~2+1)で求める軌跡はy=-x~2+1のようですがここまでのプロセスを教えてください。

  • 軌跡の描き方

    滑らかな水平面上に直線Lがあり、直線L上にOをとる。点Oを原点とし、水平面内で点Oを軸にして一定の角速度ωで上からみて反時計回りに回転するx, y軸をとる。時刻t=0では直線Lとx軸が一致しているものとする。いま、時刻t=0に点Oから距離Rだけ離れたL上の点から小球を、水平面に立っている観測者からみて、点Oに向かって 速さvoで打ち出した。xy平面内では時刻t=0の小球の位置は(x,y)=(-R,0)である。なお必要ならば小さい角δ[rad]に対して成り立つ近似式 sinδ=δ, cosδ=1を用いよ。 小球を打ち出してから微小時間Δt後の小球のy座標は y=(R-voΔt)sin(ωΔt)となっているが、ωΔtが微小であることを用いて、y≒(R-voΔt)ωΔt と近似できる。これは小球を打ち出した直後には、y方向について、等速加速度運動と近所できることを表している。そのときには小球は大きさ2mvoωの慣性力をy軸の負の向きに受けていることになる。 問 0≦t≦R/voにおける小球のxy平面内での運動の軌跡を描け。軌跡には、小球の進む向きを表す矢印をつけよ。ただし、vo>Rωとする。 小球の位置(x,y)は時刻tの関数として、 x=-(R-vot)cos(ωt) y=(R-vot)sin(ωt) と表され、小球の運動の軌跡は、半径がR-votの円の一部となるので、以下の図のようになると思います。(以下の図は、半径がRの円の一部を描き、その円周を8等分して、それぞれの点から原点Oに引いた線を8等分し、R/8づつ原点に近づいていくように点をうって、その点を結んでいく、というような手順で描きました。) しかし、解答に載っていた図は、x=-R/2に関して左右対称でした。 厳密に描くとすれば、以下の図の方が正しいと思うのですが、大学入試ではどちらでもいいのでしょうか? (これは京大実戦模試の問題です。)

  • 軌跡

    放物線y=f(x)の頂点の座標が(2a-2,-4a^2+12a-8)のときの頂点の軌跡を求める問題で x=2a-2.y=-4a^2+12a-8とおくのは x座標が2a-2.y座標が-4a^2+12a-8だからということでいいのですが 軌跡が求められるのはなぜですか? 簡単にでいいので教えてください。

  • 4次のルンゲクッタ法について

    物理か数学か,カテゴリーで悩んだのですがこちらでおねがいします. 数値積分を用いて質点の軌道予測を行いたいと思っているのですが,それには精度の高い4次のルンゲクッタを用いればよいと教えられました. 恥ずかしながらいまいちルンゲクッタを理解していないため,軌道予測ができないままでいます. 質点の軌道予測に用いるのは, 1:Δt=0.001secで採取した空力データ(風洞より) 2:採取した質点にかかる空気力は,横方向(y軸)にかかる力 3:質点はx軸方向に飛ぶ 4:F=m(d^2x/dt^2) 知りたいのは,まっすぐ飛ばした質点が,飛翔中に時々刻々と受ける横方向に力によりどれだけ曲がるか,というものです. これをルンゲクッタで解く方法を,教えていただけるととても助かります. 宜しくお願い致します.

  • 軌跡の問題が分かりません

    問題:点A(6,0)と円x~2+y~2=16上の点Qを結ぶ線分AQの中点をPとする。Qがこの円上を動くとき、点Pの軌跡を求めよ。 解:点P,Qの座標を、それぞれ(x,y),(s,t)とする。 Qは円x~2+y~2=16上にあるから s~2+y~2=16・・・(1) Pは線分AQの中点であるから x=6+s/2 y=t/2 ゆえにs=2x-6 t=2y これを(1)に代入すると (2x-6)~2+(2y)~2=16 すなわち(x-3)~2+y~2=4・・・(2) 逆に、円(2)上の任意の点は、条件を満たす。  よって、求める軌跡は、中心が(3,0)半径が2の円である。 <~2は2乗の意> 疑問:(1)にx=6+s/2 y=t/2を代入すると、なぜ点Pの軌跡が出てくるのでしょうか。よく分かりません。 よろしくお願いします。

  • 高校数学、軌跡

    tが0≦t≦1を動く時、(x-t)^2+(y-t)^2=1の軌跡を求めよ。 (x-t)^2+(y-t)^2=1は(t、t)を中心とする半径1の円である。 したがって、tが0≦t≦1を動く時、中心はy=x(0≦x≦1)を動く。ここまではわかったのですが、画像のように軌跡がなるのがわかりません。 中心が(0,0)にある半径1の円が(1,1)にある半径1の円へ平行移動したものが軌跡というのはわかるのですが、塗りつぶしている部分をどう考えているのか、黒い部分に挟まれた空白部分はなぜ生まれるのかがわかりません。教えてください。

  • 軌跡の問題で

    軌跡を求める問題で 点P(x、y)が原点を中心とする 半径1の円周上を動くとき、 点R(x(x+y)/2、y(x+y)/2)は どんな図形上を動くか という問題で 私は まず円の式はx^2+y^2=1で R(u、v)とおいて 円の式とu=x(x+y)/2、v=y(x+y)/2から 2(u+v)=x^2+y^2+2xy 2xy=2(u+v)-1・・・(1) それとは別に 2x^2=1+2(u-v)・・・(2) 2y^2=1+2(v-u)・・・(3) が分かり (2)×(3)=(1)^2から・・・・(4) uとvの関係が分かり Rの軌跡は円 2(x^2+y^2)-(x+y)=0 と言うことがわかりました しかし、答えを見ると (1)(2)(3)(4)から逆に このようなu、vについては -1≦2(u+v)≦1 となるから、(1)(2)(3)を満たすx、yの実数値が 存在する。 の一文が追加されています、この意味とどのような 求め方でこの不等式が出てきたのわかりません どなたかわかる方教えて頂けないでしょうか

  • オイラー法、ルンゲクッタ法について。

    オイラー法、ルンゲクッタ法について。 この2つについて分からない事があるので質問します。 まず、オイラーについてですが、yi+1=yi+hf(x,y)という式がテイラー展開によって求まると言われましたが、テイラー展開の2次以降の項は微少量として無視できるのは分かります。でもそもそもテイラー展開ってひとつ先の値を今の値から求まるみたいな展開でしたっけ??というのが一つ目の質問です。 2つ目は、オイラーの式の中のf(x,y)についてです。簡単なバネ・マス・ダンパ系を考えた時、運動方程式はm・d2x/dt2+c・dx/dt+kx=0となると思いますが、この場合のf(x,y)はどうやって求めるのでしょうか。 3つ目はルンゲクッタそもそもについてです。 ルンゲクッタとはK1K2K3K4という係数(?)に1221という重みをかけるとyi+1が求まるそうですが、この理由がどんなサイトや本を見ても納得出来ません。 何か分かりやすい本やサイトがあれば教えて頂けないでしょうか。 以上3つの質問、回答よろしくお願いします。

  • ルンゲクッタ法(RK4)で斜方投射の問題を解く

    下のように自由落下する物体の速度の変化を求めるプログラムをC言語で作成しました。これを応用して斜方投射された物体のx、y座標の「位置」の変化を求めるプログラムをルンゲクッタ法(RK4)で作成したいです。x,yについて2つの式を立てないといけないと思うのですがどのような式を立ててプログラムを組めかよいかよく分かりません。どなたか式とできればプログラムを教えて下さい。 k:空気抵抗、m:物体の質量 g:重力加速度としています。 #include <stdio.h> double yd(double t,double y){ double k=0.1; double m=0.1; double g=9.8; double r=g-((k*y)/m); return r; } double runge(double t,double y,double h){ double k1,k2,k3,k4,r; double h2=h/2.0; k1=yd(t,y); k2=yd(t+h2,y+(h2*k1)); k3=yd(t+h2,y+(h2*k2)); k4=yd(t+h,y+(h*k3)); r=y+(h/6.0)*((k1+(2.0*k2)+(2.0*k3)+k4)); return r; } void main(){ double y=0.0; double h=0.001; double t=0.0; int i; for(i=0;i<100;i++){ t+=h; y=runge(t,y,h); printf("%f %f\n",t,y); } }