- 微分方程式ソルバ２
微分方程式ソルバの動作がまだまだだ。ちょっとした改善をしてみた。
ndsolveとndsolve2の違い
次の例は600分割で
ndsolveで計算した例です。
setodeexpr("u1");
setodeexpr("-0.01*u1-sin(u0)");
u0[0] = 0;
u0[1] =
2.1;
nt =
600;
n_rk =
1;
t_init = 0;
t_fin =
100;
modalcolor(200,200,0);
ndsolve(2, "u0", t_init, t_fin,
nt, n_rk, "tt", "yy", [0,1], "P");
polyline("P", nt+1, "", "id");
endodeexpr();
end;
|
 | 同じ問題を分割数を減らして計算してみます。（ステップ数が減るため計算速度は一気にUPしますが結果は良くありません。）
setodeexpr("u1");
setodeexpr("-0.01*u1-sin(u0)");
u0[0] = 0;
u0[1] =
2.1;
nt =
60;
n_rk =
1;
t_init = 0;
t_fin =
100;
modalcolor(200,200,0);
ndsolve(2, "u0", t_init, t_fin,
nt, n_rk, "tt", "yy", [0,1], "P");
polyline("P", nt+1, "", "id");
endodeexpr();
end;
|
 | 良い結果が得られなかった同じ条件で
ndsolve2によって折れが少なくなるように分割ステップを調整しながら計算した結果を次に示します。ndsolve2は計算ステップを折れ線間の勾配dy/dtの差が常に
5度以下に保たれるようにステップ幅を調整しながら計算するようにしたものです。
setodeexpr("u1");
setodeexpr("-0.01*u1-sin(u0)");
u0[0] = 0;
u0[1] =
2.1;
nt =
60;
n_rk =
1;
t_init = 0;
t_fin =
100;
modalcolor(200,200,0);
nt
= ndsolve2(2, "u0", t_init, t_fin, nt, n_rk, 5, "tt", "yy", [0,1],
"P");
polyline("P", nt+1, "", "id");
endodeexpr();
end;
|
 |
|