« 2008年07月 | メイン | 2008年09月 »

2008年08月 アーカイブ

2008年08月13日

Maximaで微分方程式

せっかくMaxima使い始めたので、慣れるためにさらにMaximaのマニュアルを読み進めて遊んでると、Maximaは微分方程式も解けることがわかった。難しそうな響きを擁して実際難しい、物理学の基本ツールである。一瞬気が進まなかったが、大したものだと思ったので、敬意を表して試してみた。
すると、微分方程式が見えた気分になるようなグラフが出せることがわかり、面白かったので、ここに書き記すことにした。

まず、コマンドラインで微分させてみる。diff()が微分演算だ。

例によってTeXmacsを使って出力を整形している。(%i??)に続く青い部分が入力した部分、(%o??)に続く黒い部分がMaximaの出力だ。diff()の2つ目の引数は何で微分するか、3つ目の引数(省略可)は何回微分するかを意味する。

次に、desolve()を使ってシンプルな微分方程式\frac{d^2}{dx^2}y+\frac{d}{dx}y+y=0を解かせてみる。

desolve()の2つ目の引数は、どの関数について解くかを示す。diffの前の'(アポストロフィー)は、そのLisp関数を評価(実行&展開)せず、そのまま渡すことを指示するものである(desolve()の場合はこれを省略して評価してしまっても良いようだが、後述のode2()では不可)。y(x)の(x)はyがxの関数であることを示し、これを省略するとdy/dxが0になってしまう。
atvalue()を使って、y,y'の初期値を与えてから再度解かせてみる。

入力行の行末に;の代わりに$を付けると、計算結果を出力しない。
念のため、元の方程式に当てはまるかどうか、検算してみる。

%は直前の計算結果を表す。%o7は見るからに0なのだが、なぜか計算してくれなかったので、ratsimp()で整理させた。

常微分方程式(未知関数が1つだけの微分方程式)はode2()を使って解く方法もある。ode2()を使う場合、初期値は一般解の計算結果に対してic1(),ic2(),bc2()などを使って与える。同じ方程式を解かせてみる。

ode2()の2つ目の引数は常微分方程式の未知関数、3つ目の引数はその関数が従属する変数だ。上の場合は、引数からyがxの関数であることがわかるので、yをy(x)と書かなくても良い。

desolve()とode2()の能力を比較するため、解くには少し難しそうな微分方程式xy+\frac{d}{dx}y=0を解かせてみる。

ode2()は難なく解いた。念のため検算しておく。

それに対して、同じものをdesolve()に解かせると、下のような奇妙な式になった。

おそらく、ラプラス変換してx,yについて解いてラプラス逆変換したもの、という意味なのだろうが、これでは使えないだろう。これでも使える解なのかどうかを一応確認してみる。

やはり解けていない、使えない解であることがわかる。常微分方程式を解くのはdesolve()よりode2()の方が優れているようだ。

しかし、desolve()は複数の未知関数を含む連立微分方程式も解ける。
\left\{ \begin{array}{lll} \frac{d}{dx}f(x) & = & \frac{d}{dx}g(x)+\sin{x} \\ \frac{d^2}{dx^2}g(x) & = & \frac{d}{dx}f(x)-\cos{x} \end{array} \right.

e1とe2が方程式、f(x)とg(x)が未知関数だ。desolve()には、上のように、方程式と未知関数をそれぞれリスト形式で与える。

次に、plotdf()を使って微分方程式のグラフを出してみる。plotdf()に与える微分方程式は、dy/dx=F(x,y)または[dx/dt, dy/dt]=[F(x,y), G(x,y)]に限定されている。

load("plotdf");
plotdf(x);
これを実行すると、下のような画面が開く。(FreeBSD 6.xではTcl/Tkをインストールしていても「wishが無い」というエラーが出たが、wishという名のwish?.?へのリンクを作ればOKだった)

たくさんある矢印は、その地点(x,y)におけるdy/dxの向きと大きさを表している。このグラフ上のどこかをクリックすると、その座標を初期値とするグラフが描かれる。下の図は(0,0)辺りをクリックしたもの。

さらに色々な所をクリックする。

次の2つの図は、

plotdf(sin(x));
を実行したものである。

xとyを含む式も試してみる。

plotdf(sin(x)+cos(y));

放物線を横にした、x=y2を出そうとした。両辺をxで微分すると1=2y(dy/dx)なので、dy/dx=1/(2y)である。よりシンプルに、dy/dx=1/yで試そうとして、

plotdf(1/y);
とすると、原因がよくわからないが、Tclのエラーになった。しかし、次のように小さなオフセットを与えると、何か描画された。2つ目以降の引数に色々なオプションを指定しているが、これらが無くてもエラーにはならない。
plotdf(10^-10+1/y,[trajectory_at,0.5,1],[nsteps,200],[xradius,3],[yradius,6],[xcenter,3],[ycenter,0]);

trajectory_atで初期値を(0.5,1)として描画したものだが、左側が何かおかしい。さらに何ヶ所かクリックしてみると、

y=0の近辺でおかしくなるようだ。1/(y+⊿y)-1/yの計算でもしているのだろうか。

x=cos(t), y=sin(t)のグラフは円になるので、[dx/dt, dy/dt]=[-y,x] (=[-sin(t), cos(t)])のグラフも円になる(図は省略)。それよりもっと場を複雑にさせようと考えて、dx/dtとdy/dtを周期関数にしてみた。次のグラフは[dx/dt, dy/dt]=[sin(y),sin(x)]である。

plotdf([sin(y),sin(x)]);

色を変えながら初期値をクリックした。セル状のグラフが無限に続く。
dx/dtとdy/dtを入れ替えると、次のようなグラフになる。
plotdf([sin(x),sin(y)]);

こんなのも作ってみた。

plotdf([
10^-10 - sin(atan(y/x))*cos(%pi/4*sin(8*atan(y/x))) - cos(atan(y/x))*sin(%pi/4*sin(8*atan(y/x))),
10^-10 - sin(atan(y/x))*sin(%pi/4*sin(8*atan(y/x))) + cos(atan(y/x))*cos(%pi/4*sin(8*atan(y/x)))
]);

[dx/dt, dy/dt]=[F(x,y), G(x,y)]で表せない(dx/dtにtが入るとか)方程式も、desolve()で解けるものなら
、いつものplot2d()でグラフにできる。
\left\{\begin{array}{l}\frac{dx}{dt}=3x-7y\\ \frac{dy}{dt}=5x-4y \end{array}\right.
(この方程式は[dx/dt, dy/dt]=[F(x,y), G(x,y)]で表せるが、plotdf()ではおかしなグラフになる)

atvalue(x(t),t=0,1);
atvalue(y(t),t=0,1);
Z:desolve([diff(x(t),t)=3*x(t)-7*y(t), diff(y(t),t)=5*x(t)-4*y(t)],[x(t),y(t)]);
plot2d([parametric, x(t), y(t), [t,0,6]],[nticks,200]),Z;

アステロイド曲線\left\{\begin{array}{l} x''(t)=2y'(t)-3x(t)\\ y''(t)=-2x'(t)-3y(t) \end{array}\right.(実はx=cos(t)^3, y=sin(t)^3から逆算したもの)を描いてみる。

atvalue(x(t),t=0,1);
atvalue(y(t),t=0,0);
atvalue(diff(x(t),t),t=0,0);
atvalue(diff(y(t),t),t=0,0);
Z: desolve([
diff(x(t),t,2)=2*diff(y(t),t)-3*x(t),
diff(y(t),t,2)=-2*diff(x(t),t)-3*y(t)
],[x(t),y(t)]);
plot2d([parametric,x(t),y(t), [t,0,2*%pi]],[nticks,200]),Z;

[parametric, ...]の部分を、それらのリストにすることにより、複数のグラフを1枚に描ける。

eq: [
diff(x(t),t,2)=2*diff(y(t),t)-3*x(t),
diff(y(t),t,2)=-2*diff(x(t),t)-3*y(t)
];
atvalue(x(t), t=0, 1);
atvalue(y(t), t=0, 0);
atvalue(diff(x(t), t), t=0, 0);
atvalue(diff(y(t), t), t=0, 0);
Z1: desolve(eq, [x(t), y(t)]);
atvalue(x(t), t=0,2);
Z2: desolve(eq, [x(t), y(t)]);
atvalue(x(t), t=0,3);
Z3: desolve(eq, [x(t), y(t)]);
L1: [parametric, x(t), y(t), [t, 0, 2*%pi]], Z1;
L2: [parametric, x(t), y(t), [t, 0, 2*%pi]], Z2;
L3: [parametric, x(t), y(t), [t, 0, 2*%pi]], Z3;
plot2d([L1,L2,L3], [nticks,200]);

続きを読む "Maximaで微分方程式" »

2008年08月15日

Maximaでマンデルブロー図形

お盆なので、特にお盆とは関係ないが、ついでに定番のマンデルブロー集合も描いておく。

f(c,z):=c-z^2$
density(x,y):= block([c,z,d,i],
c:x+y*%i,
z:c,
d:0,
for i:0 while (abs(z)<2 and d<40) do (z:f(c,z),d:d+1),
d
)$
plot3d(density,[x,-0.5,1.5],[y,-1.5,1.5],[grid,100,100],
[gnuplot_pm3d,true],[gnuplot_preamble,"set view map;unset surface;set xrange [-0.5:1.5];set yrange [-1.5:1.5];"]);

続きを読む "Maximaでマンデルブロー図形" »

2008年08月17日

ペンギン in the world

Languagenotation
Japaneseペンギン
Englishpenguin
GermanPinguin
Frenchpingouin
Spanishping?ino
Italianpinguino
Simplified Chinese企鹅
Traditional Chinese企鵝
Russianпингвин
Polishpingwin
CzechTučňák
HungarianPingvin
Dutchpinguïn
Thaiเพนกวิน
Hangul펭귄
Turkishpenguen
Portuguesepinguim
Arabicالبطريق
Persianپنگوئن
Finnishpingviini
Danishpingvin
Swedishpingvin

続きを読む "ペンギン in the world" »

About 2008年08月

2008年08月にブログ「Weblog on mebius.tokaichiba.jp」に投稿されたすべてのエントリーです。過去のものから新しいものへ順番に並んでいます。

前のアーカイブは2008年07月です。

次のアーカイブは2008年09月です。

他にも多くのエントリーがあります。メインページアーカイブページも見てください。

Powered by
Movable Type 3.35