マンガでよめる痔のこと・薬のこと

オイラー法、2次ルンゲクッタ法、4次ルンゲクッタ法のC言語プログラムについて教えてください!
課題なのですが、まったくわからず困ってます><

1 常微分方程式 dy/dx=f(x,y),y(0)=1 の数値解をオイラー法を用いて計算するプログラムを作為せよ。ただし、f(x,y)=3-6x^2-4x+2xyとする。
2 α=1,β=1,γ=1/2,σ=1/2 の場合の2次ルンゲクッタ法を考える。1と同じ常微分方程式(f(x,y)も同じ)を考え、その数値解を求めるプログラムを作成せよ。また、オイラー法と2次ルンゲクッタ法の実行結果を示して、2つの近似精度を比較せよ。
3 1と同じ常微分方程式(f(x,y)も同じ)を考え、その数値解を4次ルンゲクッタ法を使って求めるプログラムを作成せよ。また、オイラー法、2次ルンゲクッタ法、4次ルンゲクッタ法の実行結果を示して、3つの近似精度を比較せよ。

以上の3つです。
休んでいた自分が悪いのですが、ネットで調べてもよくわからなくて…
わかる方、よろしくおねがいします…

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

A 回答 (2件)

下記のサイトを参照するとよいと思います。


オイラー法、2次の場合修正オイラー法などをお勉強するとよいと思います。
http://pc-physics.com/info.html
追)少しは、ネットで調べる事をお勧めします。これも勉強ですよ
    • good
    • 0

あなたのいう「オイラー法」とか「2次・4次ルンゲクッタ法」のアルゴリズムを書いてください. そして, その上でどこが分からないのかを具体的に書いてください.


「よくわからない」ということは「少しはわかる」んですよね?

この回答への補足

ネットで調べてみても、頭にハテナが浮かぶ状況ですので、「よく」ではなく「まったく」といったほうが正しいですね…
すいません…

補足日時:2010/03/04 19:19
    • good
    • 3

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

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

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

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

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

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が求められることを確かめて下さい。

Q微分方程式をルンゲクッタ法を用いて解くプログラム

2*(y''(x))+0.5(y'(x))+2.5y=1
Y(0)=0.1,y'(0)=0.1
この微分方程式をルンゲクッタ法を用いて、
x軸が0から50位(0.01刻み)までのyの値を算出するプログラムを作りたいのですが、どうすればいいのかさっぱり分かりません。
ネットで色々と調べてみたのですが、ルンゲクッタ法をどうやってプログラムで表せば良いのかが分かりません。
どうか、回答をよろしくお願いします。

Aベストアンサー

もし学校の課題なら自分で考えないと自分の将来のためによくないよ。
そうでないなら、基本はここで十分でしょう。
http://www.damp.tottori-u.ac.jp/%7Eooshida/edu/ode/
サンプルソフトが欲しければ、ここ。
http://www.unige.ch/math/folks/hairer/software.html

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ルンゲクッタ法の二階微分方程式(Fortran)

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

-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となってしまうんですがどうなんでしょう?

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

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

-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...続きを読む

Aベストアンサー

卒論で微分方程式の数値計算をやったので多少の覚えはあります。結局僕の場合は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で一発で出ますよ。グラフも。

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 を計算す...続きを読む

Q2階微分方程式の解法

d^2y/dx^2+2*x*dy/dx=0
境界条件 x=0: y=1、x→∞: y→0
上記の微分方程式をルンゲクッタ法を使って数値的に解きたいのですが、どのように x→∞ : y→0 の境界条件をいれればいいかわかりません。どなたか教えて頂けないでしょうか。回答よろしくお願いします。

Aベストアンサー

一部誤りがありましたので訂正します。
解は減少関数の筈なので次の初期値で計算を始めてみて下さい。
y=1, y'=-1
x が大きくなった場合、yの値は0 に近づく筈ですね。

y=1, y'=-0.5; y=1, y'=-2;
等でも xが∞に近づく時に y=0 の解に近づきそうなので一意に解が得られるのか、
或いは一意に解が求まる為には、初期値や制約条件が不足しているのかも知れません。
何か考え落しがありそうに思えてきました。 
以上のような考えを参考に、数値計算を続けて正しい答が得られるように頑張ってみてください。

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があるのにsqrtが・・・

最初にあるプログラムを作っていたのですが、その祭sqrtでエラーが出てしまったので、別にsqrtを使う短いプログラムを作りました。それが以下のものになります。(test3.c)

#include<stdio.h>
#include<math.h>

int main(void){
int atai;
double x=2.0;
atai=(int)sqrt(x);
printf("atai is %c",atai);
}


sqrtが動作するか確かめる為のものです。このプログラム(test3.c)で以下のエラーが出てしまいました。

/tmp/ccqEejZ1.o(.text+0x4d): In function `main':
test3.c: undefined reference to `sqrt'
collect2: ld はステータス 1 で終了しました

ちなみにLinux(Fedora core 4)を使用しています。
初心者ということもありなぜエラーが出るのか分かりません。
ちなみに<math.h>を使用しないプログラムは普通に動作します。
解決方法をご存知の方がいらっしゃいましたらご教授下さいませ。
よろしくお願いします。

最初にあるプログラムを作っていたのですが、その祭sqrtでエラーが出てしまったので、別にsqrtを使う短いプログラムを作りました。それが以下のものになります。(test3.c)

#include<stdio.h>
#include<math.h>

int main(void){
int atai;
double x=2.0;
atai=(int)sqrt(x);
printf("atai is %c",atai);
}


sqrtが動作するか確かめる為のものです。このプログラム(test3.c)で以下のエラーが出てしまいました。

/tmp/ccqEejZ1.o(.text+0x4d): In function `main':
test3.c: undefined r...続きを読む

Aベストアンサー

コンパイルのオプションつけてないとかないですよね?
たとえばTerminalとかだとコンパイル時に 『-lm』とかつけますが

Q角度の計算でdegとradがあると思いますが、これの違いやまた計算で角

角度の計算でdegとradがあると思いますが、これの違いやまた計算で角度を求める際の手順がよくわかりません。お分かりになるかたお願いします。

Aベストアンサー

小学校の算数で、分度器を使って教えられる角度の単位がdeg(degreeの略)です。つまり、直角=90°=90degです。
高校で三角関数の微積分を習うと出てくる角度単位がrad(弧度、radianの略)です。半径と弧の長さが同じ扇形の中心角を1radとした角度単位です。従って、直角=90°=π/2radです。
角度単位をradにすると、三角関数の微積分が簡単になるという長所があります。
双方の単位間の換算式は回答が既出ですから省略します。

Q差分法とオイラー法の違いについて

最近微分方程式の数値解析について学びだした者です。

微分方程式の数値解法として差分法とオイラー法があると思うのですが、この2つの違いや互いの位置づけはどうなっているのでしょうか?

また、差分法には風上法などがありますが、これらとオイラー法の位置づけについても教えていただきたいです。

できればこれら近辺の全体的な体系について教えていただけるとうれしいです。よろしくお願いします。

Aベストアンサー

 すいません、#1です。#1では話を簡単にしようと思って、差分法の定式化を簡略にし過ぎました。もう気づいていると思いますが、

 オイラー法:
  y[i+1]=(1+h)*y[i]     (1)

に対して、差分法:
  y[i+1]-(1+h)*y[i]=0   (2)

としてしまっては、(1)と(2)は全く同じです。差分法で、

 (y[i+1]-y[i])/h=(y[i+1]+y[i])/2   (3)

くらいの事はしておかないと、領域型の性質を持ちません。領域型をめざすなら、これ以上は簡単にできないと思います。(3)と(1),(2)の違いですが、(1),(2)では(移行すれば)、

 (y[i+1]-y[i])/h=y[i]    (4)

で、(3)との違いは、右辺にy[i+1]を考慮するかどうかだけです。

 (3)では左辺の数値微分の結果を、x[i+1]とx[i]の間の中点における微分係数と解釈し、それを中点における関数値に等しいとしています。中点における関数値は、y[i+1]とy[i]の平均で、十分良く近似できるとも仮定しています。

 (1),(2),(4)では、y[i+1]を計算するために、y[i]の情報しか使いませんが、(3)では(今は一点だけですが)まわりの点の情報も使って、平均的に良好な解を得ようとします。ここが直接法と領域型の発想の違いだと言えます。

 式の上ではほんのわずかな違いですが、結果の違いはけっこうあります。添付図は、Excelで計算したものです(^^;)。

  ・黒線が、y=e^x
  ・青丸が、(3)
  ・赤線が、(1),(2),(4)

の結果です。

 すいません、#1です。#1では話を簡単にしようと思って、差分法の定式化を簡略にし過ぎました。もう気づいていると思いますが、

 オイラー法:
  y[i+1]=(1+h)*y[i]     (1)

に対して、差分法:
  y[i+1]-(1+h)*y[i]=0   (2)

としてしまっては、(1)と(2)は全く同じです。差分法で、

 (y[i+1]-y[i])/h=(y[i+1]+y[i])/2   (3)

くらいの事はしておかないと、領域型の性質を持ちません。領域型をめざすなら、これ以上は簡単にできないと思います。(3)と(1),(2)の違いですが、(1),(2)...続きを読む


人気Q&Aランキング