Euler法
をテンプレートにして作成
home
>
サイトマップ
開始行:
RIGHT:寄稿:東條遼平
* 常微分方程式とは [#u9491d2d]
常微分方程式とは一変数関数の導関数が含まれている方程式の...
例えば次の式であれば、yの導関数がtとyによって記述されてい...
#ref(img100.png,nolink)
多くの物理現象は微分方程式によって記述され、
例を挙げれば、放射性物質の崩壊、空気抵抗、電気回路、減衰...
そのため微分方程式を解く事は大変重要なのですが、簡単な微...
数式で表示された解である解析解を求めるのは大変困難になり...
そこでPCを使い数値解を求めてみようと思います。
* Euler法とは [#i1b3ea9b]
#ref(img101.png,nolink)
このような方程式を考えてみます。y'(t)とはy(t)の傾きの事で...
初期値y(0)、y'(0)が判っているのであれば、
十分小さい幅dtを考える事でy(dt)がわかり、y(dt)とdtを使って
y'(dt)が求められます。それを繰り返す事で任意のtでのy(t)を...
イメージとしては次のようになります。
#ref(euler1.png,nolink)
これを式で書くと、
#ref(img102.png,nolink)
であり、t=tnのとき、y(t)=yn、y'(t)=y'n
とすると、yn+1は次のように表せます。
#ref(img103.png,nolink)
これを使って微分方程式の数値解を求めるのがEuler法です。
* 連立微分方程式 [#dc813872]
連立微分方程式とはその名の通り、ある関数の導関数に他の関...
例えば次のようなものがあります。
#ref(img107.png,nolink)
ここで、
#ref(img108.png,nolink)
のようにおくと、
#ref(img109.png,nolink)
と出来ることから、連立微分方程式も引数と変数が増えるだけ...
わかります。
* 高次微分方程式 [#fc86b40e]
#ref(img110.png,nolink)
のように高次の導関数を含む微分方程式を考えてみます。この...
#ref(img111.png,nolink)
というふうに式変形をし、
#ref(img112.png,nolink)
とすると、これは上で書いた連立微分方程式と同じ形になり、y...
初期値がわかれば同じ方法で解けることになります。
つまり、
#ref(img100.png,nolink)
という一番シンプルな微分方程式さえ解ければ、連立微分方程...
高次微分方程式も解けるということです。
* Taylor展開における一次近似 [#i7bcdf7b]
また、無限回微分出来るという事から、Taylor展開を考えると、
#ref(img104.png,nolink)
#ref(img105.png,nolink)
#ref(img106.png,nolink)
となり、2階微分以降を無視するとEuler法の式と同じになりま...
また、dt2以降を削る事になるため、誤差はdtの二乗のオーダー...
精度はかなり悪いと言えます。
* ソースコード [#cd285df3]
素直にEuler法のソースコードを書くならば、tとy(t)の値を引...
定義して、その関数ポインタfとt0からtnまでの区間をdiv分割...
次のようになります。
double Euler(double (*f)(double t, double y),
double t0, double y, double tn, int div)
{
double h, t;
int i;
h = (tn - t0);
if(div) h /= div;
t = t0;
for(i=0; i<div; i++){
y += h*(*f)(t, y);
t += h;
}
return y;
}
これでy(tn)の値を返してくれるのですが、これでは
二階の微分方程式や、連立微分方程式を解こうとすると
新たに受け取る引数と初期値の個数が増えただけで同じ処理を...
作らなければいけません。そこで関数ポインタと初期値の数を...
汎用な関数を作る事にします。
void Euler(double (*f[])(double t, double *x),
double t0, double *x, double tn, int ...
{
double h, t;
double *temp;
int i, j;
if((temp = (double*)malloc(sizeof(double)*num)) == NUL...
fprintf(stderr, "Allocation Error\n");
exit(1);
}
h = (tn - t0);
if(div) h /= div;
t = t0;
for(i=0; i<div; i++){
for(j=0; j<num; j++){
temp[j] = h*(*f[j])(t, x);
}
for(j=0; j<num; j++){
x[j] += temp[j];
}
t += h;
}
free(temp);
}
関数ポインタの配列と初期値を配列で渡し、その要素数をnumで...
連立微分方程式の変数がいくつになっても、また、何次の導関...
この関数で解を求めることができます。例えばf(t,x,y,z)であ...
x[0]をxに、x[1]をyに、x[2]をzに対応させます。
double x[] = {1, 0};
double (*f[2])(double, double*);
f[0] = f1;
f[1] = f2;
Euler(f, 0, x, 10, 1000, 2);
printf("%f %f\n", x[0], x[1]);
のようにすればt=10のときの値だけを求められますし、
int i;
double t=0, tn=10;
double x[] = {1, 0};
double (*f[2])(double, double*);
double h;
f[0] = f1;
f[1] = f2;
h=(tn-t)/1000;
printf("%f %f\n", t, x[0]);
for(i=0; i<1000; i++){
Euler(f, t, x, t+h, 1, 2);
t += h;
printf("%f %f\n", x[0], x[1]);
}
のようにすれば、グラフ作成用のデータも取ることができて便...
- &ref(main.c);
- &ref(calculation.c);
- &ref(calculation.h);
終了行:
RIGHT:寄稿:東條遼平
* 常微分方程式とは [#u9491d2d]
常微分方程式とは一変数関数の導関数が含まれている方程式の...
例えば次の式であれば、yの導関数がtとyによって記述されてい...
#ref(img100.png,nolink)
多くの物理現象は微分方程式によって記述され、
例を挙げれば、放射性物質の崩壊、空気抵抗、電気回路、減衰...
そのため微分方程式を解く事は大変重要なのですが、簡単な微...
数式で表示された解である解析解を求めるのは大変困難になり...
そこでPCを使い数値解を求めてみようと思います。
* Euler法とは [#i1b3ea9b]
#ref(img101.png,nolink)
このような方程式を考えてみます。y'(t)とはy(t)の傾きの事で...
初期値y(0)、y'(0)が判っているのであれば、
十分小さい幅dtを考える事でy(dt)がわかり、y(dt)とdtを使って
y'(dt)が求められます。それを繰り返す事で任意のtでのy(t)を...
イメージとしては次のようになります。
#ref(euler1.png,nolink)
これを式で書くと、
#ref(img102.png,nolink)
であり、t=tnのとき、y(t)=yn、y'(t)=y'n
とすると、yn+1は次のように表せます。
#ref(img103.png,nolink)
これを使って微分方程式の数値解を求めるのがEuler法です。
* 連立微分方程式 [#dc813872]
連立微分方程式とはその名の通り、ある関数の導関数に他の関...
例えば次のようなものがあります。
#ref(img107.png,nolink)
ここで、
#ref(img108.png,nolink)
のようにおくと、
#ref(img109.png,nolink)
と出来ることから、連立微分方程式も引数と変数が増えるだけ...
わかります。
* 高次微分方程式 [#fc86b40e]
#ref(img110.png,nolink)
のように高次の導関数を含む微分方程式を考えてみます。この...
#ref(img111.png,nolink)
というふうに式変形をし、
#ref(img112.png,nolink)
とすると、これは上で書いた連立微分方程式と同じ形になり、y...
初期値がわかれば同じ方法で解けることになります。
つまり、
#ref(img100.png,nolink)
という一番シンプルな微分方程式さえ解ければ、連立微分方程...
高次微分方程式も解けるということです。
* Taylor展開における一次近似 [#i7bcdf7b]
また、無限回微分出来るという事から、Taylor展開を考えると、
#ref(img104.png,nolink)
#ref(img105.png,nolink)
#ref(img106.png,nolink)
となり、2階微分以降を無視するとEuler法の式と同じになりま...
また、dt2以降を削る事になるため、誤差はdtの二乗のオーダー...
精度はかなり悪いと言えます。
* ソースコード [#cd285df3]
素直にEuler法のソースコードを書くならば、tとy(t)の値を引...
定義して、その関数ポインタfとt0からtnまでの区間をdiv分割...
次のようになります。
double Euler(double (*f)(double t, double y),
double t0, double y, double tn, int div)
{
double h, t;
int i;
h = (tn - t0);
if(div) h /= div;
t = t0;
for(i=0; i<div; i++){
y += h*(*f)(t, y);
t += h;
}
return y;
}
これでy(tn)の値を返してくれるのですが、これでは
二階の微分方程式や、連立微分方程式を解こうとすると
新たに受け取る引数と初期値の個数が増えただけで同じ処理を...
作らなければいけません。そこで関数ポインタと初期値の数を...
汎用な関数を作る事にします。
void Euler(double (*f[])(double t, double *x),
double t0, double *x, double tn, int ...
{
double h, t;
double *temp;
int i, j;
if((temp = (double*)malloc(sizeof(double)*num)) == NUL...
fprintf(stderr, "Allocation Error\n");
exit(1);
}
h = (tn - t0);
if(div) h /= div;
t = t0;
for(i=0; i<div; i++){
for(j=0; j<num; j++){
temp[j] = h*(*f[j])(t, x);
}
for(j=0; j<num; j++){
x[j] += temp[j];
}
t += h;
}
free(temp);
}
関数ポインタの配列と初期値を配列で渡し、その要素数をnumで...
連立微分方程式の変数がいくつになっても、また、何次の導関...
この関数で解を求めることができます。例えばf(t,x,y,z)であ...
x[0]をxに、x[1]をyに、x[2]をzに対応させます。
double x[] = {1, 0};
double (*f[2])(double, double*);
f[0] = f1;
f[1] = f2;
Euler(f, 0, x, 10, 1000, 2);
printf("%f %f\n", x[0], x[1]);
のようにすればt=10のときの値だけを求められますし、
int i;
double t=0, tn=10;
double x[] = {1, 0};
double (*f[2])(double, double*);
double h;
f[0] = f1;
f[1] = f2;
h=(tn-t)/1000;
printf("%f %f\n", t, x[0]);
for(i=0; i<1000; i++){
Euler(f, t, x, t+h, 1, 2);
t += h;
printf("%f %f\n", x[0], x[1]);
}
のようにすれば、グラフ作成用のデータも取ることができて便...
- &ref(main.c);
- &ref(calculation.c);
- &ref(calculation.h);
ページ名:
home
>
Modified by
物理のかぎプロジェクト
PukiWiki 1.4.5_1
Copyright © 2001-2005
PukiWiki Developers Team
. License is
GPL
.
Based on "PukiWiki" 1.3 by
yu-ji
Powered by PHP 5.3.29HTML convert time to 0.002 sec.