とらりもん - C言語による微分方程式の解析 (Euler法) Diff
- Added parts are displayed like this.
- Deleted parts are displayed
like this.
dx/dt = 2x をEuler法で解こう(t =0のときx =1とする)。課題6-1と同じ解析である:
/* euler1.c
* 2004/02/15 K. Nishida
*/
# include <stdio.h>
main()
{
double t=0;
double dt=0.1;
double x=1;
double dx;
printf("%lf %lf\n", t, x);
for (t=dt; t<=1;t=t+dt)
{dx=2*x*dt;
t=t+dt){
dx=2*x*dt;
x=x+dx;
printf("%lf %lf\n", t, x);
}
}
}
注意: 左端の数字(6.など)は打ち込まないで良い。
$ gcc euler1.c -o euler1
$ ./euler1 > euler1.txt
$ gnuplot
GNUPLOT> plot "euler1.txt" w p, exp(2*x) w l
{{attach_view(euler1.png)}}
上の図で, 緑色の線が解析解(x=exp(2t)), 赤点がオイラー法による数値解である。
'''課題13-1''' 上を参考にして, dx/dt = 2x (1- x) をEuler法で解け(t =0のときx =0.1とする)。課題5-2と同じ解析である。
'''課題13-2''' 以下のLotka-Volterra方程式を解き, 結果をGNUPLOTで表示せよ:
dS/dt = pS - qWS
dW/dt = -rW + sWS
t =0でS =2, W =0.5とする。p =q =r =s =1とする。
!運動方程式
Newtonの運動方程式は,
d2x / dt2 = F/m (xは質点の座標, Fは力, mは質量)
であり, 2階の微分方程式なので, 直接的には積分できない。そこで, v=dx/dtという変数を導入し, xとvの2つについての1階の連立微分方程式に変換する:
dv / dt = F/m
dx / dt = v
あとは, 上のLotka-Volterra方程式などと同様にEuler法で解くことができる。
たとえば, バネの振動運動を解いてみよう。バネ定数kのバネに質量mのおもりがついて振動しているときの運動方程式は,
dv / dt = - x * k/m
dx / dt = v
となる。初期条件および, パラメータ(ここではmとk)を適当に設定して解いてみよう。以下にその例を示す。
/* euler2.c
* spring motion
* 2004/02/15 K. Nishida
*/
# include <stdio.h>
main()
{main(){
double t=0;
double dt=0.01;
double x=1, v=0;
double dx, dv;
double k=1;
double m=1;
printf("%lf %lf %lf\n", t, x, v);
for (t=dt; t<=10; t=t+dt)
{dx=v*dt;
dv=-1*x*(k/m)*dt;
x=x+dx;
v=v+dv;
printf("%lf %lf %lf\n", t, x, v);
for (t=dt; t<=10; t=t+dt){
dx=v*dt;
dv=-1*x*(k/m)*dt;
x=x+dx;
v=v+dv;
printf("%lf %lf %lf\n", t, x, v);
}
}
注意: 左端の数字(6.など)は打ち込まないで良い。
$ gcc euler2.c -o euler2
$ ./euler2 > euler2.txt
$ gnuplot
GNUPLOT> plot "euler2.txt" using 1:2 title "x" w l, "euler2.txt" using 1:3 title "v" w l
{{attach_view(11b.png)}}
'''課題13-3''' おもりに空気抵抗が働く場合を考えて, 上で作ったバネの振動のプログラムを改造せよ。
ヒント:空気抵抗は,
dv / dt = - x * k/m - v * a (aは抵抗係数というパラメータ)
という形で運動方程式に入る。aに適当な値を設定して18行目をそのように改造すれば良い。
以下にこの課題の解の例を示す:
{{attach_view(11c.png)}}
'''課題13-4''' 前問を解析的に解いて, その解析解と, 前問の数値解をいっしょにグラフに表示してくらべてみよ。
/* euler1.c
* 2004/02/15 K. Nishida
*/
# include <stdio.h>
main()
{
double t=0;
double dt=0.1;
double x=1;
double dx;
printf("%lf %lf\n", t, x);
for (t=dt; t<=1;
{dx=2*x*dt;
dx=2*x*dt;
x=x+dx;
printf("%lf %lf\n", t, x);
}
}
注意: 左端の数字(6.など)は打ち込まないで良い。
$ gcc euler1.c -o euler1
$ ./euler1 > euler1.txt
$ gnuplot
GNUPLOT> plot "euler1.txt" w p, exp(2*x) w l
{{attach_view(euler1.png)}}
上の図で, 緑色の線が解析解(x=exp(2t)), 赤点がオイラー法による数値解である。
'''課題13-1''' 上を参考にして, dx/dt = 2x (1- x) をEuler法で解け(t =0のときx =0.1とする)。課題5-2と同じ解析である。
'''課題13-2''' 以下のLotka-Volterra方程式を解き, 結果をGNUPLOTで表示せよ:
dS/dt = pS - qWS
dW/dt = -rW + sWS
t =0でS =2, W =0.5とする。p =q =r =s =1とする。
!運動方程式
Newtonの運動方程式は,
d2x / dt2 = F/m (xは質点の座標, Fは力, mは質量)
であり, 2階の微分方程式なので, 直接的には積分できない。そこで, v=dx/dtという変数を導入し, xとvの2つについての1階の連立微分方程式に変換する:
dv / dt = F/m
あとは, 上のLotka-Volterra方程式などと同様にEuler法で解くことができる。
たとえば, バネの振動運動を解いてみよう。バネ定数kのバネに質量mのおもりがついて振動しているときの運動方程式は,
dv / dt = - x * k/m
となる。初期条件および, パラメータ(ここではmとk)を適当に設定して解いてみよう。以下にその例を示す。
/* euler2.c
* spring motion
* 2004/02/15 K. Nishida
*/
# include <stdio.h>
{
double t=0;
double dt=0.01;
double x=1, v=0;
double dx, dv;
double k=1;
double m=1;
for (t=dt; t<=10; t=t+dt)
{dx=v*dt;
dv=-1*x*(k/m)*dt;
x=x+dx;
v=v+dv;
for (t=dt; t<=10; t=t+dt){
dx=v*dt;
dv=-1*x*(k/m)*dt;
x=x+dx;
v=v+dv;
printf("%lf %lf %lf\n", t, x, v);
}
}
注意: 左端の数字(6.など)は打ち込まないで良い。
$ gcc euler2.c -o euler2
$ ./euler2 > euler2.txt
$ gnuplot
GNUPLOT> plot "euler2.txt" using 1:2 title "x" w l, "euler2.txt" using 1:3 title "v" w l
{{attach_view(11b.png)}}
'''課題13-3''' おもりに空気抵抗が働く場合を考えて, 上で作ったバネの振動のプログラムを改造せよ。
ヒント:空気抵抗は,
dv / dt = - x * k/m - v * a (aは抵抗係数というパラメータ)
という形で運動方程式に入る。aに適当な値を設定して18行目をそのように改造すれば良い。
以下にこの課題の解の例を示す:
{{attach_view(11c.png)}}
'''課題13-4''' 前問を解析的に解いて, その解析解と, 前問の数値解をいっしょにグラフに表示してくらべてみよ。