数値計算の演習問題で以下の二階微分方程式をルンゲ・クッタで解けという問題があります。

-y"+x^2・y=e・y(eは定数、”・”は単なる掛け算)

y(0)=1, y'(0)=0, 0<=x<=2までを計算せよ。



これは

y’=z・・・(1)
 
z'=(x^2-e)y・・・(2)

この2つの連立方程式を解けばよいところまではわかります。
まず(2)を解くときにルンゲ・クッタの場合
(k1+2k2+2k3+k4)/6の項(←公式の右辺第二項)のk(1~4)を求めなければいけません。

質問はkの求め方です。

本にはy'=f(x,y,z) , z'=g(x,y,z)とおけば
(2)の場合だと例えばk1は
k1=g(xn,yn,zn)dxで計算する。と書いてあります
しかしz'=(x^2-e)y(←zが入ってない) なので、計算すると
k1=g(xn,yn)dxとなってしまうんですがどうなんでしょう?

おそらくどこかで勘違いしてると思うんです。
長い質問になってしまいましたがどうかご教授のほどよろしくお願いします。

A 回答 (2件)

卒論で微分方程式の数値計算をやったので多少の覚えはあります。

結局僕の場合はy'の項が入ってるのでリープフロッグの方が収束が早いという事でそっちにしましたが。

> k1=g(xn,yn)dxとなってしまうんですがどうなんでしょう?

問題ありません。

> 本にはy'=f(x,y,z) , z'=g(x,y,z)とおけば

というのはy', z'がx,y,zの関数という一般の場合であり、仮に
f(x, y, z) = x + y
という関数が存在しても何の問題もないのと同じです。

k1を計算する際にも
k1 = (x^2 - e)y * dx
と値が一意に決まり、問題ないでしょう。

ちなみにこの程度の微分方程式なら刻み幅dxにも寄りますが、プログラムなど組まなくてもExcelで一発で出ますよ。グラフも。
    • good
    • 0
この回答へのお礼

どうもありがとうございます。
もしかしたら私は勘違いしてるんじゃないかと思ってたので・・
しつこく質問してしまいすいませんでした。

お礼日時:2001/07/07 03:53

良いと思いますよ。


ちなみに、私は、dx=hとして、
k1=f(x,y,z)・h
l1=g(x,y,z)・h
k2=f(x+h/2,y+k1/2,z+l1/2)・h
l2=g(x+h/2,y+k1/2,z+l1/2)・h
k3=f(x+h/2,y+k2/2,z+l2/2)・h
l3=g(x+h/2,y+k2/2,z+l2/2)・h
k4=f(x+h,y+k3,z+l3)・h
l4=g(x+h,y+k3,z+l3)・h
としてプログラムを書きました。
    • good
    • 0
この回答へのお礼

どうもありがとうございます。
同じ問題やられてたんですか(笑)
お答えの中で
k1=f(x,y,z)・h
l1=g(x,y,z)・h
とありますが
実際プログラムで計算するとき
f(x,y,z)の中身は
f(x,y,z)=z

g(x,y,z)の中身は
g(x,y,z)=(x^2-e)・y
ということでいいんですか?

って何がいいたいかと言うと
例えばg(x,y,z)・hを実際に計算する場合、
g(x,y,z)=(x^2-e)・y
に変数(xとy)に数を代入して計算しますよね。
その場合g(x,y,z)(←zが入っている)にはどこにもzを代入する
ところがないんですがg(x,y,z)と記述するのは単に形式的なものと
考えていいんですか?

お礼日時:2001/07/06 00:30

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

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

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

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

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エクセルで計算すると2.43E-19などと表示される。Eとは何ですか?

よろしくお願いします。
エクセルの回帰分析をすると有意水準で2.43E-19などと表示されますが
Eとは何でしょうか?

また、回帰分析の数字の意味が良く分からないのですが、
皆さんは独学されましたか?それとも講座などをうけたのでしょうか?

回帰分析でR2(決定係数)しかみていないのですが
どうすれば回帰分析が分かるようになるのでしょうか?
本を読んだのですがいまいち難しくて分かりません。
教えてください。
よろしくお願いします。

Aベストアンサー

★回答
・最初に『回帰分析』をここで説明するのは少し大変なので『E』のみ説明します。
・回答者 No.1 ~ No.3 さんと同じく『指数表記』の『Exponent』ですよ。
・『指数』って分かりますか?
・10→1.0E+1(1.0×10の1乗)→×10倍
・100→1.0E+2(1.0×10の2乗)→×100倍
・1000→1.0E+3(1.0×10の3乗)→×1000倍
・0.1→1.0E-1(1.0×1/10の1乗)→×1/10倍→÷10
・0.01→1.0E-2(1.0×1/10の2乗)→×1/100倍→÷100
・0.001→1.0E-3(1.0×1/10の3乗)→×1/1000倍→÷1000
・になります。ようするに 10 を n 乗すると元の数字になるための指数表記のことですよ。
・よって、『2.43E-19』とは?
 2.43×1/(10の19乗)で、
 2.43×1/10000000000000000000となり、
 2.43×0.0000000000000000001だから、
 0.000000000000000000243という数値を意味します。

補足:
・E+数値は 10、100、1000 という大きい数を表します。
・E-数値は 0.1、0.01、0.001 という小さい数を表します。
・数学では『2.43×10』の次に、小さい数字で上に『19』と表示します。→http://ja.wikipedia.org/wiki/%E6%8C%87%E6%95%B0%E8%A1%A8%E8%A8%98
・最後に『回帰分析』とは何?下の『参考URL』をどうぞ。→『数学』カテゴリで質問してみては?

参考URL:http://ja.wikipedia.org/wiki/%E5%9B%9E%E5%B8%B0%E5%88%86%E6%9E%90

★回答
・最初に『回帰分析』をここで説明するのは少し大変なので『E』のみ説明します。
・回答者 No.1 ~ No.3 さんと同じく『指数表記』の『Exponent』ですよ。
・『指数』って分かりますか?
・10→1.0E+1(1.0×10の1乗)→×10倍
・100→1.0E+2(1.0×10の2乗)→×100倍
・1000→1.0E+3(1.0×10の3乗)→×1000倍
・0.1→1.0E-1(1.0×1/10の1乗)→×1/10倍→÷10
・0.01→1.0E-2(1.0×1/10の2乗)→×1/100倍→÷100
・0.001→1.0E-3(1.0×1/10の3乗)→×1/1000倍→÷1000
・になります。ようするに 10 を n 乗すると元の数字になるた...続きを読む

Q線形・非線形って何ですか?

既に同じようなテーマで質問が出ておりますが、
再度お聞きしたく質問します。

※既に出ている質問
『質問:線形、非線型ってどういう意味ですか?』
http://oshiete1.goo.ne.jp/kotaeru.php3?q=285400
結局これを読んでもいまいちピンと来なかった...(--;


1.線形と非線形について教えてください。
2.何の為にそのような考え方(分け方)をするのか教えてください。


勝手なお願いですが、以下の点に留意いただけると大変うれしいです。
何せ数学はそんなに得意ではない人間+歳なので...(~~;

・わかりやすく教えてください。(小学生に説明するつもりぐらいだとありがたいです)
・例をあげてください。(こちらも小学生でもわかるような例をいただけると助かります)
・数式はなるべく少なくしてください。

『そんな条件じゃ説明できないよー』という方もいると思いますが、どうぞよろしくお願いいたしますm(__)m

Aベストアンサー

昨日「線形の方がなんとなくてわかりやすくないですか」と書いたんですが、やっぱり理系の人間らしく、もうちょっときちんと説明してみます。昨日は数式をなるべく出さないように説明しようとがんばったんですが、今日は少しだけ出しますが、勘弁してください。m(__)m(あと、長文も勘弁してください)


数学的にはちょっとここまで言えるかわかりませんが、自然界の法則としては、「線形」が重要な意味を持つのは、xの値が変化するにつれて変化するyがあったときに、

(yの増加量)/(xの増加量)=A(一定)

という規則が成り立つからです。

xやyの例としては昨日の例で言う例1だとxがガムの個数、yが全体の金額、例2だとxが時間、yが走った距離です。

この規則が何で役に立つかというと、式をちょっと変形すると、

(yの増加量)=A×(xの増加量)・・(1)

ということがわかります。つまり、Aの値さえわかれば、xが増えたときのyの値が容易に推測できるようになるわけです。


ここで「Aの値さえわかれば」と書いていますが、この意味を今から説明します。

自然界の法則を調べるためには何らかの実験を行います。例えば、りんごが木から落ちる運動の測定を行います。
ここから質問者様がイメージできるかわかりませんが、りんごは時間が経つにつれて(下に落ちるにつれて)落下するスピードが速くなるんです。今、実験として、1秒ごとにりんごのスピードを測定したとします。そしてその結果をグラフにプロットしていくと、直線になることがわかります。(ここがわかりにくいかもしれませんが、実際に実験を行うとそのようになるのです)

数学の問題のように初めから「時速100kmで走る」とか「1個100円のガム」とかいうことが与えられていれば直線になることはすぐにわかります。
しかし、自然界の法則はそうもうまくいきません。つまり、実験を行ってその結果をプロットした結果が直線状になっていたときに初めて「何らかの法則があるのではないか」ということがわかり、上で書いた「Aの値さえわかれば」の「A」の値がプロットが直線状になった結果、初めてわかるのです。

そして、プロットが直線状になっているということは、永遠にそうなることが予想されます。つまり、今現在はりんごが木から落ちたときしか実験できませんが、その結果を用いて、もしりんごが雲の上から落としたときに地面ではどのくらいのスピードになるかが推測できるようになるわけです。ここで、このことがなぜ推測できるようになるかというと、(1)で書いた関係式があるからです。このように「なんらかの法則があることが推測でき、それを用いて別の事象が予言できるようになる」ことが「線形」が重要だと考えられる理由です。

しかし、実際に飛行機に乗って雲の上からりんごを落としたらここで推測した値にはならないのです。スカイダイビングを想像するとわかると思いますが、最初はどんどんスピードが上がっていきますが、ある程度でスピードは変わらなくなります。(ずっとスピードが増え続けたら、たぶんあんなに空中で動く余裕はないでしょうか??)つまり、「線形から外れる」のです。

では、なぜスピードが変わらなくなるかというと、お分かりになると思いますが、空気抵抗があるからなんですね。(これが昨日「世の中そううまくはいかない」と書いた理由です)つまり、初めは「線形」かと思われたりんごを落とすという実験は実際には「非線形」なんです。非線形のときは(1)の関係式が成り立たないので、線形のときほど容易には現象の予測ができないことがわかると思います。


では、非線形だと、全てのことにおいて現象の予測が難しいのでしょうか?実はそうでもありません。例えば、logは非線形だということをNo.5さんが書かれていますが、「片対数グラフ」というちょっと特殊な形のグラフを用いるとlogや指数関数のグラフも直線になるんです。つまり、普通のグラフでプロットしたときに「非線形」になるため一見何の法則もないように見えがちな実験結果が「片対数グラフ」を用いると、プロット結果が「線形」になってlogや指数関数の性質を持つことが容易にわかり、それを用いて現象の予測を行うことが(もちろん単なる線形よりは難しいですが)できるようになるわけです。


これが私の「線形」「非線形」の理解です。つまり、

1) 線形の結果の場合は同様の他の事象の推測が容易
2) 非線形の場合は同様の他の事象の推測が困難
3) しかし、一見非線形に見えるものも特殊な見方をすると線形になることがあり、その場合は事象の推測が容易である

このことからいろいろな実験結果は「なるべく線形にならないか」ということを目標に頑張ります。しかし、実際には先ほどの空気抵抗の例のように、どうしても線形にはならない事象の方が世の中多いんです。(つまり、非線形のものが多いんです)

わかりやすいかどうかよくわかりませんが、これが「線形」「非線形」を分ける理由だと思っています。

やっぱり、「線形の方がなんとなくわかりやすい」くらいの理解の方がよかったですかね(^^;;

昨日「線形の方がなんとなくてわかりやすくないですか」と書いたんですが、やっぱり理系の人間らしく、もうちょっときちんと説明してみます。昨日は数式をなるべく出さないように説明しようとがんばったんですが、今日は少しだけ出しますが、勘弁してください。m(__)m(あと、長文も勘弁してください)


数学的にはちょっとここまで言えるかわかりませんが、自然界の法則としては、「線形」が重要な意味を持つのは、xの値が変化するにつれて変化するyがあったときに、

(yの増加量)/(xの増加量)=...続きを読む

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

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

Aベストアンサー

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

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

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

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

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

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

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

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

頑張って下さい.

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

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

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

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

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

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

Qfortranでのcsvファイルを出力する方法

fortranでコンマ切りテキストファイル(CSVファイル)を
出力する方法をご存知の方、いらっしゃいませんでしょうか?


私が使用した処理系(intel fortran 8.0 for Windows)で

write(600,*)1.4,',',2.33,',',30.11

というコードを実行させると
fort.600 というファイル名のテキストファイルが生成されて
中身は

___1.400000_____,___2.330000_____,___30.11000____

となります。( _ は半角空白を意味します。)
私が欲しいのは余計な空白の入っていない

1.4,2.33,3.11

または

1.400000,2.330000,30.11000

という出力なのですが、どのような書式指定をすれば
これが可能なのでしょうか?

全くのfortran初心者です。よろしくご教示ください。

Aベストアンサー

汎用的な方法は、#2さんのように文字列として書き出して、adjustlやtrimを使って空白を取るのでしょうが、Intel Fortran8.0ならば以下のようなformatでできます。
write(600,10) 1.4,',',2.33,',',30.11
10 format(f0.6,a1,f0.6,a1,f0.6)
マイナス記号があっても問題ないです。書式幅が0の指定はIとFに使えます。なかなか便利ですよ。ちなみにConpaq Visual Fortranでも使えます。

Qセグメンテーション違反

C言語を使用しています。

構造体に値をいれようとしたら、コンパイルは出来るのですが、実行時に
「セグメンテーション違反です (core dumped)」
となってしまい、それ以上行えません。

構造体と代入したい変数との型は、合っています。

いろいろ本などで見ましたが、何が原因かわからず困っています。
教えてください。
宜しくお願いします。

Aベストアンサー

OSは何でしょうか。コンパイラは何を使用していますか?
通常、デバッグオプションをつけて実行すると、異常の発生したソースの箇所で止まりますので、それが手がかりになります。またNo1の方が言われてますように、ソースが公開できるのであれば、ソースを提示するのが良いかと思います。

Q高階連立常微分方程式の数値計算

4次のルンゲクッタ法を用いた数値計算を勉強しています.
1階連立常微分方程式と高階常微分方程式は理解でき,プログラムも作成することができました.

次に高階の連立常微分方程式を解こうと思ったら,頭が混乱してしまいました.

4次のルンゲクッタ法を用いて高解連立常微分方程式を解く考え方を教えて頂ければ嬉しいです.
また何か良い参考書があれば教えて頂きたいと思います.

よろしくお願いします.

Aベストアンサー

すべての変数の再高階でない微係数を、あらたな変数にしてしまえば、どんな高階の連立常微分方程式も1階連立常微分方程式になります。

たとえば、
x''+2y''+3x'y'=0
x''+4x''+5x'+6y'+7x+8y=0
とかだったら(超適当です)、
z=x' , w=y'
とすれば、
z'+2w'+3zw=0
z'+4w'+5z+6w+7x+8y=0
z-x'=0
w-y'=0
と4変数の1階連立常微分方程式になります。

QfortranでのNaNについて

お忙しい所、NaNについて回答お願い致します。

fortranでプログラムを使用して数値計算を行っていますが、計算途中で結果がNaNとなり困っています。
NaNは、0で除算を行なったり、負の数の指数を求めようとした時に出る無限大や定義できない値のようなのですが、式を見ても0で除算を行なったり、負の数の指数を求めようとしている所はありません。

もし同じような経験をされて見事解決された方、居られましたらその時の対処法をお聞かせ願えないでしょうか。
何卒よろしくお願い致します。

Aベストアンサー

環境(OSとコンパイラ)がわからないのでできるとは断言できないのですが、
浮動小数点レジスタの制御フラグをいじって、NaNが生成されたときに
例外が起きるようにしておけば、どこでNaNになるのかが絞れるのではないでしょうか?

いずれにしろ、注意深く値の変化を追っていくよりないと思います。

参考URL:http://www.xlsoft.com/jp/products/intel/cvf/docs/vf-html/pg/pg22.htm

Qファイル出力の場所を指定

現在C++にてhtmlファイルを出力するプログラムを作っているのですが、出力場所を指定することはできるのでしょうか?(現在はそのプログラムソースが保存されている場所と同じファイル内に出力されますが、それをデスクトップに出力するなど。)
もし、方法がありましたら、教えてください。
ソースや参考HPのURLなどのせていただけたらありがたいです。
環境はVisualStudio.NET2003です。
よろしくお願いします。

Aベストアンサー

単にファイル名の前にパスを指定する。

絶対パス指定
fp=fopen("c:/temp/test.txt","w");

相対パス指定
fp=fopen("./hoge/test.txt","w");


デスクトップはOSやユーザによって場所が異なるので、少し面倒です。
XPの場合環境変数を利用してこんな感じで出来ると思います。

例:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

void main(void)
{
FILE *fp;
char fname[1024];
strcpy(fname,getenv("USERPROFILE"));
strcat(fname,"/デスクトップ/test.txt");
fp=fopen(fname,"w");
//処理
fclose(fp);
}

QFORTRAN 複素数 (実部を取り出す方法)

fortran77で複素数の実部を取り出したいのですが、
関数が存在するかわかりません。どなたか、ご存知の方、
教えてください。虚部は出せました。
FORTRAN77(SALFORD FTN77)を使っています。

確か、fortran90なら real(○+□i)だったと思います。

Aベストアンサー

real() でいいと思いますが。

complex x/(1,2)/
write(*,*) x
write(*,*) real(x)
write(*,*) imag(x)
end


このQ&Aを見た人がよく見るQ&A

人気Q&Aランキング