【先着1,000名様!】1,000円分をプレゼント!

コンピュータ(プログラミング)のカテゴリに投稿しようかとも考えましたが、物理学に関する数値計算ということなので、物理学のカテゴリに投稿しました。

単振動や減衰振動、そして強制振動などを数値解析したいと思い、ルンゲクッタ法を使いシミュレートしてみようとしています。ルンゲクッタ法という方法を全く知らなかったので、インターネットや図書館で調べたのですが、どうしても分からないことろがあるので質問することにしました。

書籍やネットの情報を参考にしながら、単振動の場合を数値解析してみました(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'で表されるので、この式を変形(?)したものが、入るのかな、という予想程度にしか分かりません。
ソースをどのように変えれば減衰振動を解析できるのでしょうか。

どなたか詳しい方、教えてください。お願いします。

このQ&Aに関連する最新のQ&A

A 回答 (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)
の形にするだけ.
    • good
    • 2

・f1 や f2 は何を表す関数ですか?


・その引数である t, x, v はそれぞれ何を意味するのですか?

この回答への補足

すみませんでした。
重要なことを書き忘れていたようです。

f1,f2,t,x,vはそれぞれ、
f1:dx/dt
f2:dv/dt
t:時間
x:変位
v:速度
を表しています。
単振動の運動方程式において、出てくる、質量mとバネ定数kはどちらも1と置きました。

補足日時:2008/07/11 16:43
    • good
    • 1

このQ&Aに関連する人気のQ&A

お探しのQ&Aが見つからない時は、教えて!gooで質問しましょう!

このQ&Aを見た人はこんなQ&Aも見ています

このQ&Aを見た人が検索しているワード

このQ&Aと関連する良く見られている質問

Q振り子をオイラー法で

プログラムで、振り子のヒモの長さ、重さ、角度を自由に設定し、振り子のアニメーションが表示されるようにしようと考えています。
普通の運動方程式での振り子の運動の表示はできたのですが、
オイラー法やルンゲクッタ法を使って導き出した振り子の動きも表示したいと考えています。
ところが、どうやってオイラー法を振り子の運動方程式に応用すれば、振り子のX、y座標が分かるのか?ということが、勉強不足でどうにもわかりません…

どなたか、ヒントだけでも教えていただけると幸いです。

Aベストアンサー

 #1です。
 補足を拝見しました。

 運動方程式で、途中で速度をかませてよいのでしたら、オイラー法や修正オイラー法を使うことができます。この方法は、常に、速度⇒変位の順に解いていくものです。
http://www.asahi-net.or.jp/~hy9n-knk/eq_simu.htm

 そのほかに、邪道かもしれませんが、速度の計算を一度だけで省略する方法もなくはありません。
 2階微分方程式を差分方程式の形に直すと、初期変位として2点が必要になります。1点は問題自体に与えられているものを使いますが、2点目は、速度をかませて求めた変位を使います。あとは常に最初の2点が確保されますので、順に数値計算を進めていくことができます。

 なお、オイラー法を使う場合は、誤差がかなり大きくなってくるので注意してください。
http://maya.phys.kyushu-u.ac.jp/~knomura/education/numerical-physics/text5/node2.html

 また、ルンゲ・クッタ法でも解くことができます。
 運動方程式に適用した例を掲載したサイトがありましたので紹介します。
http://www6.ocn.ne.jp/~simuphys/runge-kutta.html

 #1です。
 補足を拝見しました。

 運動方程式で、途中で速度をかませてよいのでしたら、オイラー法や修正オイラー法を使うことができます。この方法は、常に、速度⇒変位の順に解いていくものです。
http://www.asahi-net.or.jp/~hy9n-knk/eq_simu.htm

 そのほかに、邪道かもしれませんが、速度の計算を一度だけで省略する方法もなくはありません。
 2階微分方程式を差分方程式の形に直すと、初期変位として2点が必要になります。1点は問題自体に与えられているものを使いますが、2点目は、速度を...続きを読む

Q三角関数の記述の仕方

タイトルそのまんまなんですが、三角関数はC言語ではどのように記述すればいいでしょうか?
角度にラジアン表記でπ(パイ)を使いたいんですが、その表記方法もわかりません。
僕の持っている本に載ってなかったので質問させていただきました。
よろしくお願いします。

Aベストアンサー

C言語で三角関数を使うためには、math.h をインクルードする必要があります。使い方は例えば、こんな感じです。

#define M_PI 3.14159265358979 /* 円周率 */

double x, y, theta;

theta = M_PI / 4.0;
x = cos(theta); /* sin,cos,tanの引数は弧度法の角度です。*/
y = sin(theta);

πは上記の例のように自分で定義して使ってください。

Q方程式を2分法を用いて解くプログラム

学校で出されたCプログラムの課題で、1問だけどうしても出来ない問題があるんです。
「方程式 f(x) = x2 - 2 = 0 を 2 分法を用いて解くプログラムを作成せよ。ここで、方程式 f(x) = x2 - 2 は関数として定義せよ。上位の方から 4 桁目まで正しい値が出たらループを止めるようにする。」
というものなのですが、この「2分法」というやり方もよく分かりません。
プログラムの作成方法と併せて教えて頂けると幸いです。

Aベストアンサー

こんな感じでしょうか?

#include <stdio.h>

/* 関数 f(x) */
double f(double x) {
 return x*x-2.0;
}

/* 二分法 初期値 x1<x2 と 誤差限界 eps を入力 */
double bisec(double x1, double x2, double eps) {
 double x;
 while (x2 - x1 >= eps) {
  x = (x1+x2)/2.0; /* 中点計算 */
  if (f(x1)*f(x) > 0.0) { /* 同符号か判定 */
   x1 = x;
  } else {
   x2 = x;
  }
 }
 return (x1+x2)/2.0;
}

int main(void) {
 double eps=0.00001;
 printf("%lf %lf\n",bisec(-2,0,eps), bisec(0,2,eps));

 return 0;
}

Q物理シミュレーションをする時、どのくらいプログラミングの知識があればいいのか?

物理学で自然現象をパソコンでシミュレーションするとき、プログラム言語はどのようなものを使うのでしょうか?よく知られたC言語やJava等は使わないのでしょうか?専門的なプログラム言語がいろいろあるのでしょうか?

今後パソコンで物理シミュレーションを行うことになったとき、プログラムに関してはどのくらいの知識があればいいのでしょうか?
基本的な本を見るとBasicやFortranを使ってシミュレーションの説明をしているものがありますが、このような基本的な言語も使えるようになったほうがよいのでしょうか?
C言語やJava、VBなど一般的によく知られたプログラミング言語も覚えたほうがよいのでしょうか?

Aベストアンサー

いくつかのケースがあります.

どなたかもすでにご回答のとおり,
「シミュレーション記述言語」というようなものが
いくつかありまして,
mathematica
freefem
その他,有料のfemソフトとか
流体解析ソフト

に入力するデータを作って入れる
ということになります.
それを一種の言語を用いたプログラミングだといえばそうでしょう.

別のケースとして,Cなどのプログラム言語で,
femソフトを作る場合もあります.
そのときに,どういう手順をたどるかというと,
物理現象を,マトリックス形式にする式展開の載っている本を購入します.
その本は,ハンドブックだったりしますが,
そこに載っているマトリックス形式の式を
プログラムします.
そのときには,プログラム言語は何でもいいと思います.

一部で,Fortranや,BASICというのがありますが,
これらは,Fortranについて言えば,昔ながらのものであり,これを用いたプログラムがたくさんあるという理由です.BASICも同様です.

もしも新たに,プログラムについて覚えるのならば,
CかC++でしょうね.

これからは,プログラマーは,適当に書いて,
最適化コンパイラが汗をしぼって「高速で効率的な実行形式を生成する」時代です.

私も,言語はBASICから入り,大学時代はFortran
を主に使ってました.シミュレーションをプログラムするレベルなら,どれも大差ありません.
Cが一番小回りが利いていいと思いますね.

いくつかのケースがあります.

どなたかもすでにご回答のとおり,
「シミュレーション記述言語」というようなものが
いくつかありまして,
mathematica
freefem
その他,有料のfemソフトとか
流体解析ソフト

に入力するデータを作って入れる
ということになります.
それを一種の言語を用いたプログラミングだといえばそうでしょう.

別のケースとして,Cなどのプログラム言語で,
femソフトを作る場合もあります.
そのときに,どういう手順をたどるかというと,
物理現象を,マトリックス形...続きを読む

Q単振り子の運動方程式

重力加速度g、質量m、紐の長さl、空気抵抗無視。

単振り子の運動方程式はこうなりますよね。
mlθ"=-mgsinθ
これがよくわからないのです。
どういう座標系についての運動方程式なのですか?

軌道にそってx軸を定めると
θl=x
mx"=-mgsinθ  軌道に沿った運動方程式?
⇔mlθ"=-mgsinθ  どういう座標系の運動方程式なの?
そしてこれの一般解はどういう風になりますか?
初期条件としてt=0でθ=φとします。

Aベストアンサー

まず座標系についてのお話をします。下の図をご覧下さい。

  y
  ↑
  ・→x
   \
   →\
   θ \
      ●

振子の支点を・、先端に吊るされたおもりを●で表しています。支点の位置をxy座標の原点に取るならば、鉛直からの振れ角をθとして
x= l sinθ  (1)
y= -l cosθ  (2)
であることは既にご承知かと思います。
このように置くこと自体が、(x, y)の直交座標系から(l, θ)の極座標系に移行していることに相当します。ただほとんど自明なことなので「極座標に置き換えて」などとわざわざ断っていないわけです。
極座標系に移行したことで問題の本質はx(t), y(t)の代わりにl(t), θ(t)を求めることに帰着します。大抵の場合はひもは伸び縮みしないと仮定しますのでlについて解く必要はなく、θについてのみ解くことになります。その方程式が
ml(d^2θ/dt^2)= -mg sinθ  (3)
なわけです。

しかしこの方程式は初等関数の範囲では解くことが出来ません。そこで初等物理の範囲ではθが小さい場合に限って問題を考えることにし、
sinθ≒θ  (4)
の近似を行って解きます。このとき(3)は
ml(d^2θ/dt^2) = -mg θ  (5)
となります。これの解き方はいろいろあります。線形微分方程式の理論を知っていれば解は直ちに
θ= C sin{√(g/l) t+α} ←Cは定数  (6)
だと分かります。αはC sinα=φを満たす定数です。
2階の微分方程式ですが初期条件が「t=0でθ=φ」の一つしか与えられていないので、定数が一つ未定のまま残ります(*1)。

愚直に微分方程式を解くのであれば下のようにやります。
l(d^2θ/dt^2)(dθ/dt) = -g θ(dθ/dt)
d/dt {(dθ/dt)^2} = -(g/l) d/dt (θ^2) ←両辺に(dθ/dt)をかけた上で、積の導関数の公式((y^2)'=2y y')を逆に使った
(dθ/dt)^2 = -(g/l) θ^2 +C1 ←C1は積分定数
dθ/dt = √{-(g/l) θ^2 +C1}  (7)
ここでθ=√(l/g)√C1 sinψと変数を変換すると
dθ/dt = √C1√(1-sin^2 ψ)  (8)
を経て
√(l/g)√C1 cosψ dψ = √C1 cosψ dt  (9)
と変形でき、両辺を積分することで
√(l/g) ψ= t+C2 ←C2は積分定数  (10)
を得ます。θの表式に戻すと
θ=√(l/g)√C1 sin{√(l/g) (t+C2)}  (11)
となります。これは本質的に(6)と同じ式です。初期条件「t=0でθ=φ」を代入することで
φ=√(l/g)√C1 sin{√(l/g)C2}  (12)
を得ます。これを使うと(11)からC1, C2のいずれかを消去できます。初期条件がもう一つあれば運動は一意に定まります(脚注参照)。

もちろん、「軌道に沿ってx軸を定める」でも解けます。この場合の運動方程式は
m(d^2 x/dt^2)= -mg sin(x/l)  (13)
となります。本質的に(3)と同じであることは申し上げるまでもなく、同様に解くことができます。

考え方は上記でよいはずですが中間で計算ミスがあるかも知れませんので、ONEONEさんご自身でも確認しながら読んで頂けると幸いです。

*1 もし初期条件が「t=0でθ=φまでおもりを持ち上げて手を放す」という意味であれば、「θの最大値はφ(厳密には|φ|)」という条件が新たに加わるので運動は一意に定まります。この場合はφsinα=φからα=π/2、よってθ=φsin{√(g/l) t+(π/2)}=φcos{√(g/l) t}と求めることができます。

まず座標系についてのお話をします。下の図をご覧下さい。

  y
  ↑
  ・→x
   \
   →\
   θ \
      ●

振子の支点を・、先端に吊るされたおもりを●で表しています。支点の位置をxy座標の原点に取るならば、鉛直からの振れ角をθとして
x= l sinθ  (1)
y= -l cosθ  (2)
であることは既にご承知かと思います。
このように置くこと自体が、(x, y)の直交座標系から(l, θ)の極座標系に移行していることに相当します。ただほとんど自明なことなので「極座標に置き換えて」...続きを読む

Qコンパイルエラー invalid operands to binary

自己啓発で入力文字列をBASE64デコードする関数を作っているのですが、L20~L23(a[0] = strchr(b64, p[0]) - b64;)でコンパイルエラーinvalid operands to binaryが発生して色々試行錯誤しているのですが、どうしてもエラーがとれません。
ソースをここに書くのは大変恐縮なのですが、原因がわかる方がいらっしゃいましたら、教えていただけないでしょうか?

char *Base64n(unsigned char *buf, size_t length, size_t *outlen)
{
const char b64[] = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrst         uvwxyz0123456789+/=";
unsigned char *p;
unsigned char *q;
unsigned char a[4];
char *RtnBuf;
int j=0;
int cnt;

RtnBuf = (char *)malloc(length+1);
memset(RtnBuf, 0, length+1);

p = (unsigned char*)buf;
q = (unsigned char*)RtnBuf;
cnt = 0;
while(*p != 0) {
a[0] = a[1] = a[2] = a[3] = 0;
a[0] = strchr(b64, p[0]) - b64;
a[1] = strchr(b64, p[1]) - b64;
a[2] = strchr(b64, p[2]) - b64;
a[3] = strchr(b64, p[3]) - b64;

q[0] = ((a[0] << 2) | (a[1] >> 4)) & 0xff;
cnt++;

if (p[2] != '=') {
q[1] = ((a[1] << 4) | (a[2] >> 2)) &0xff;
cnt++;
}
if (p[3] != '=') {
q[2] = ((a[2] << 6) | a[3]) & 0xff;
cnt++;
}
p += 4;
q += 3;
}
*outlen = cnt;
return(RtnBuf);
}

コンパイルはRed Hatでgccを使ってコンパイルしています。
引数は第1引数がデコード対象の文字列、第2引数がデコード対象文字列長、第3引数がデコード後の文字列長で、戻り値がデコード後の文字列です。

自己啓発で入力文字列をBASE64デコードする関数を作っているのですが、L20~L23(a[0] = strchr(b64, p[0]) - b64;)でコンパイルエラーinvalid operands to binaryが発生して色々試行錯誤しているのですが、どうしてもエラーがとれません。
ソースをここに書くのは大変恐縮なのですが、原因がわかる方がいらっしゃいましたら、教えていただけないでしょうか?

char *Base64n(unsigned char *buf, size_t length, size_t *outlen)
{
const char b64[] = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrst ...続きを読む

Aベストアンサー

手元にgccが無いので、Windowsのbcc32でコンパイルしてみたところ、
エラーが出ずに通ってしまいました。
だから自信がないのですが…。

a[0] = strchr(b64, p[0]) - b64;



a[0] = strchr(b64, p[0]) - b64[0];

に変えてみたらどうでしょう。
"invalid operands to binary"
は、たぶん、マイナスの両側で型が違っていることを
表しているのではないでしょうか。(←これも自信なし)
strchr()が返すのはchar *型ですが、
b64の型は、charの配列型です。
型が違うので、コンパイラが不正と判断したのかもしれません。

とするとbcc32でなぜ通ったかが問題になるのですが…。
C言語規格でも、ポインタ同士の引き算のところは
ややこしくなっています。
規格解釈の違いがあるのかもしれません。

Q【C言語】二階微分方程式をルンゲクッタで解く解き方が…

二階微分方程式をルンゲクッタで解くプログラムを作っているのですが、上手く合成関数が定義できず、上手く行きません。途中までプログラムを作ったので、見ていただけませんか??

問題
http://12lien.web.fc2.com/q.jpg

(PC環境Windows でX on windows 3 コンパイラはgcc )
 
プログラム⇒http://12lien.web.fc2.com/base.txt

プログラム中ではdx/dt = Y, x = X とおきました。

オネガイシマスm(_ _)m

Aベストアンサー

えっと, 2階微分方程式 d^2x /dt^2 = f(dx/dt, x) を Runge-Kutta で数値的に解くために dx/dt = y とおいて連立 1階微分方程式d(y, x)/dt = (f(y, x), y) を解きにいってるんですよね?
少なくとも, プログラム上で dy/dt を計算する関数 FUN に x を渡さないとダメですし, dx/dt を計算する関数は不要 (単に y を返すだけだし, それなら関数を呼ぶ必要はない) です.
でどうするかなんですが, まず dy/dt を計算する関数
double FUN(double t, double y, double x) {
ここで t, x, y を使って dy/dt を計算する
}
を作ります. そして, 実際に Runge-Kutta で解く関数では 1ステップずつ「これを呼出して y[i] を逐次計算⇒その情報を基にx[i] を計算」を繰り返すことになります. 挙げられたプログラムでは 2重ループになってるけど, これは無駄です. 大雑把には
for (i = 0; i < 1000; ++i) {
FUN を 4回呼出して t, x[i-1], y[i-1] から y[i] と x[i] を計算
}
という感じ. なお, FUN の 2回目以降の呼出しのときに x[i] もそれらしい値を使うよう注意する必要があります.

えっと, 2階微分方程式 d^2x /dt^2 = f(dx/dt, x) を Runge-Kutta で数値的に解くために dx/dt = y とおいて連立 1階微分方程式d(y, x)/dt = (f(y, x), y) を解きにいってるんですよね?
少なくとも, プログラム上で dy/dt を計算する関数 FUN に x を渡さないとダメですし, dx/dt を計算する関数は不要 (単に y を返すだけだし, それなら関数を呼ぶ必要はない) です.
でどうするかなんですが, まず dy/dt を計算する関数
double FUN(double t, double y, double x) {
ここで t, x, y を使って dy/dt を計算す...続きを読む

Qニュートン法を使って解を求めるC言語プログラム

C言語を使って y=x^2-4x のyの解をニュートン法を使って求める
プログラムを作る課題を出されたんですが、ニュートン法が良く分かっていないので、いろいろ調べたり、人に聞いたりしたところ
#include<stdio.h>
#include<math.h>
void main()
{
int counter=0;
double an,g,f,sh=0.0001;
printf("初期値を入力して下さい==>");
scanf("%ld",&an);
do{
g=(an*an)/(2*an-4);
f=2*an-4;
counter++;
}while(fabs(f)>sh);
printf("反復回数 %d 回 y=%lf \n",counter,g);
}
でプログラムがこんな感じになったんですが、結局ニュートン法がどうなのかがわかりません。 なんか微分とかやるとか言われたんですが、工業系の学校で数学の授業が無いので微分についてがわかりません。
このプログラムは、コンパイルはできるんですが、動きません。
ニュートン法についてよくわからないのでどこが間違ってるかわかりません。 ニュートン法についてできるだけ分かりやすく解説してほしいです。

C言語を使って y=x^2-4x のyの解をニュートン法を使って求める
プログラムを作る課題を出されたんですが、ニュートン法が良く分かっていないので、いろいろ調べたり、人に聞いたりしたところ
#include<stdio.h>
#include<math.h>
void main()
{
int counter=0;
double an,g,f,sh=0.0001;
printf("初期値を入力して下さい==>");
scanf("%ld",&an);
do{
g=(an*an)/(2*an-4);
f=2*an-4;
counter++;
}while(fabs(f)>sh);
printf("反復回数 %d 回 y=%lf \n",counter,g);
}
でプログラムがこんな感...続きを読む

Aベストアンサー

ニュートン法を理解するには、せめて微分の初歩の初歩は理解して下さい。
とりあえず、必要なところだけ。

y=x^2-4x
の解をもとめるには、次のようにします。
計算途中の近似解をa(n)とすると、それを元に得られる、より精度の良い近似解a(n+1)は、
a(n+1)=(a(n)*a(n))/(2*a(n)-4)
として求められます。(ニュートン法の理論から)
これを繰り返していって、
|a(n+1)-a(n)| < 0.0001
となったときに、十分精度の良い解が得られたと判断し、計算を終了します。
(計算終了の閾値0.0001は提示されたプログラムから取りました。)

プログラムの間違いは、下記の2点。
誤:scanf("%ld",&an);
正:scanf("%lf",&an);

誤:f=2*an-4;
正:f=g-an; an=g;

上記を修正し、初期値が2より大きい場合は4.000000が、初期値が2未満のときは
0.000000が求められることを確かめて下さい。

Qe^-2xの積分

e^-2xの積分はどうしたらよいのでしょうか…。e^xやe^2xsinxなどはのってるのですがこれが見つかりません。お願いします。

Aベストアンサー

いささか、思い違いのようです。

e^-2x は、 t=-2x と置いて置換してもよいけれど、牛刀の感がします。

e^-2x を微分すると、(-2)*( e^-2x )となるので、

e^-2x の積分は、(-1/2)*( e^-2x )と判明します。

Qオイラー法とルンゲ・クッタ法

「オイラー法とルンゲ・クッタ法の計算精度を数値的に比較しなさい」と課題を出されましたがさっぱりわからないのです.
最低でもオイラー法と2次のルンゲ・クッタ法を比較しないといけないのですがどのような方法でどのような結果になるのでしょうか?
お願いします

Aベストアンサー

「数値的に」ということは,理論的な解析はやらなくても良いようですね.

直接的にやり方をお教えするのは禁止されていますので,簡単に手順の概要をしますから,あとは実際にやって考えてみて下さい.

1. 数値計算に頼らなくても,解が解析的に求まる微分方程式を用意する
(調和振動子,バネとかが良いかもしれません)

2. 「刻み幅」を同じ値にし,オイラー法とルンゲクッタ法で解を求めていく

3. 解析解(微分方程式を実際に解いた式)と数値計算結果の差を比較する
(グラフを描いてみると良いかもしれません.差の絶対値を取る方が良いと思います)

4. 「刻み幅」(ステップサイズ)の大きさを半分,もしくは2倍にして,上記1~3を再度試みる

5. 時間の許す限り,4. を繰り返す
(刻み幅を10倍とかやってみても良いかもしれません)

参考URLの本に,この辺のことは全部書いていますので,御購入なさってはいかがでしょうか.

頑張って下さい.

参考URL:http://tinyurl.com/65udeo

「数値的に」ということは,理論的な解析はやらなくても良いようですね.

直接的にやり方をお教えするのは禁止されていますので,簡単に手順の概要をしますから,あとは実際にやって考えてみて下さい.

1. 数値計算に頼らなくても,解が解析的に求まる微分方程式を用意する
(調和振動子,バネとかが良いかもしれません)

2. 「刻み幅」を同じ値にし,オイラー法とルンゲクッタ法で解を求めていく

3. 解析解(微分方程式を実際に解いた式)と数値計算結果の差を比較する
(グラフを描いてみる...続きを読む


人気Q&Aランキング