![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
chdir("path") - ワーキングディレクトリの変更
pwd - ワーキングディレクトリの表示
exit - scilabを終了する quit - 1段階上のレベルへ移動する。トップレベルだとscilabを終了する。
%i 虚数単位 %pi 円周率 %e 自然対数の底 %inf 無限大 %eps 機械イプシロン %nan %inf-%inf
%s = poly(0,'s') %z = poly(0,'z') $ = poly(0,'$')
%t - 真 %f - 偽
predefコマンドを使って $SCIHOME/.scilabに定数を定義することができる。
Vine-4.1のscilabではSCIHOMEは $HOME/.Scilab/scilab-3.1.1/ と定義されている。
+ 和 - 差 * 積(行列同士の場合は行列の積) / 商(行列同士の場合は A/B=A*inv(B)となる。) ^ べき(行列の場合は行列のべき) \ 左からの割り算(スカラーの場合はA\B=(1/A)*B、行列の場合はA\B=inv(A)*B)
.* 成分同士の積 ./ 成分同士の商 .\ 成分同士の左からの割り算 .^ 成分個々のべき
== 等しい < より小さい > より大きい <= 以上 >= 以下 <> or ~= 等しくない
& 論理積 | 論理和 ~ 論理否定
スラッシュ2つ(//)以降はコメントと見なされる。
コマンドの区切りは改行かコンマかセミコロン。
セミコロンだと結果は表示されない。
長いコマンドを次の行に続けて書く場合は行の終わりにピリオド３つ(...)をつける。
if 式 then
コマンド;
:
elseif 式 then
コマンド;
：
else
コマンド;
:
end
thenの代わりに改行かコンマでもよい。
select 式
case 値 then
コマンド;
:
case 値 then
コマンド;
:
else
コマンド;
:
end
thenの代わりに改行かコンマでもよい。
for 制御変数 = 範囲式 do
コマンド;
コマンド;
:
end
doのかわりに改行かコンマでもよい。 範囲式は 1:10 とか 1:0.1:1 など。範囲式がリストの場合は 各要素についてループを実行し、行列の場合は各列ベクトルについて ループを実行する。
while 条件式 do
コマンド
コマンド;
:
end
doのかわりに改行かコンマでも良い。
for文やwhile文から抜け出すにはbreakが使える。
continueはループ内のそれ以後の処理をスキップし、次のループまで飛ぶ。
input(message) - キーボードから数値や変数名を入力 input(message,"string") - キーボードから文字列を入力 disp( val ) - 任意のオブジェクトを決められたフォーマットで画面に出力 printf(format, args) - C言語と同様のprintf関数を使って画面に出力 mprintf(format,matrix) - matrixの各行毎にprintf関数を使って画面に出力
fp=file('open',filename) - ファイルオープン(obsolate)
file('close',fp) - ファイルクローズ(obsolate)
fscanf(fp,format,val,val,...) - ファイル読み込み(obsolate)
fprintf(fp,format,val,val,...) - ファイル書き込み(obsolate)
fd=mopen(filename,"wt") - ファイルオープン ( w : 書き込み, r : 読み込み, t : テキスト, b : バイナリ ) mclose(fd) - ファイルクローズ mfprintf(fd,format,val,val,....) - ファイル書き込み mfprintf(fd,format,matrix) - matrixの各行毎にformatに従ってファイル書き込み mfscanf(fd,format) - 戻り値はベクトル
fscanfMat(filename) - 数値行列を一気に読み込む fprintfMat(filename,mat) - 数値行列を一気に書き込む
その他、mput,mputl,mputstr,mget,mgetl,mgetstrなど。
◎ファイルから行列を読み込む例
M=fscanfMat("mat0.txt")
mat0.txtの中身 -----------
#
# comment
#
1 2 3
4 5 6
7 8 9
--------------------------
◎ファイルへ行列を書き込む例
fprintfMat("mat1.txt",M [,format] )
formatはC言語のprintfの書式(デフォルトは"%f")
d=dir('.'); カレントディレクトリ取得
d.name ファイル名のリスト
d.date 更新時間のリスト
d.isdir ファイルがディレクトリかどうかのフラグのリスト
listfiles(path) : ファイル名のリストを返す。
プログラムをファイルに記述し、それを呼び出して実行したい場合は
exec( "path" [,mode] )
とする。プログラムファイルの拡張子は".sce"とすることが多い。
modeオプション
0 : デフォルト値 -1 : 何も表示されない 1 : 各コマンドラインエコー 2 : プロンプト --> が表示される 3 : エコー + プロンプト 4 : 各プロンプト前に中断．キャリッジリターン後に実行再開 7 : 中断 + プロンプト + エコー : デモに役立つモード
関数のみを記録した外部ファイルから関数を読み込むには
getf("path" [,opt])
とする。関数ファイルの拡張子は".sci"とすることが多い。 optは
"c" : 関数をコンパイルする。(デフォルト) "n" : 関数をコンパイルしない。 "p" : 関数をプロファイリング用にコンパイルする。
注) 最近は関数の読み込みもexecを使うことが推奨されているらしい。
これらの基本関数は行列に対しては要素毎の計算となる。
基本演算
abs - 絶対値, 大きさ ceil - 小数点切上げ conj - 複素共役 cumprod - 累積積 cumsum - 累積和 round - 丸め込み fix - ゼロの方向に丸める floor - 切り捨て frexp - 浮動小数点を２つの成分、指数と仮数に分解する int - 整数部分 real - 実数部 imag - 虚数部 imult - 虚数単位 i の掛け算 log - 自然対数 log10 - 対数 log1p - 1を加えられた引数の自然対数の計算 log2 - 底が２の対数 minus - (-) 減算演算子、符号変換 modulo - m を法とする余り pmodulo - m を法とする正の余り size - オブジェクトのサイズ sign - 符号関数 sqrt - 平方根
三角関数
sin - 正弦 sinh - 双曲線正弦 cos - 余弦関数 cosh - 双曲線余弦関数 acos - 要素毎の逆余弦 acosh - (要素毎の)逆双曲線余弦 asin - 逆正弦 asinh - 逆双曲線正弦 atan - 第2象限と第4象限の逆正接 atanh - 逆双曲線正接 cotg - 余接 (cotangent) coth - 双曲線余接 tan - 正接 tanh - 双曲線正接
その他の関数
amell - ヤコビアン楕円積分関数 besseli - 第一種修正ベッセル関数 (I sub alpha). besselj - 第一種ベッセル関数 (J sub alpha). besselk - 第二種修正ベッセル関数 (K sub alpha). bessely - 第二種ベッセル関数 (Y sub alpha). beta - ベータ関数 gamma - ガンマ関数 gammaln - ガンマ関数の対数 binomial - 二項分布確率 erf - 誤差関数 erfc - 相補誤差関数 erfcx - スケーリングされた相補誤差関数 calerf - 誤差関数の計算 dlgamma - gammaln 関数、ψ 関数の微分 legendre - 準ルジャンドル関数
直接配列を生成するには
[1,2,3,4]等 - 行(row)ベクトル [1;2;3;4]等 - 列(column)ベクトル
範囲とステップを指定して生成するには
x1:dx:x2 - x1からx2までをdxステップで刻んだ配列(行ベクトル)
linspace(x1,x2,n) - x1からx2をn点に等分割した配列(行ベクトル)を生成。
0から10まで1間隔ならlinspace(0,10,11)とする。
logspace(x1,x2,n) - 対数的に等間隔の配列(行ベクトル)を生成。
生成した行ベクトルを列ベクトルに変換するには(')を使って転置する。
ベクトルは行数あるいは列数が1の行列と見なされる。
直接行列を生成するには
[1,2,3; 3,4,5; 6,7,8]
特殊な行列
zeros(m,n) - m行n列のゼロ行列 ones(m,n) - 全部1の行列(m行n列) eye(n,n) - nxnの単位行列 rand(m,n,'uniform') - [0,1)の乱数行列 rand(m,n,'norm') - 正規分布乱数行列
○行列Aと同じサイズの行列を作る場合eye(A)やrand(A)などとすることができる。
○複数の行列から C=[A,B]などのようにして新しい行列を作ることができる。(縦に並べるときはC=[A;B]とする)
(入出力の項参照)
コロン(:)で範囲指定ができる。コロンのみだと行、あるいは列全体を示す。
インデックスは1から始まる。
A(m,n) - m,n要素 A(m,:) - m行 A(:,n) - n列 A(m1:m2,n) - n列のm1行からm2行まで。列の指定も同様 A(m1:m2,n1:n2) - m1,m2,n1,n2で指定される範囲 A($) - 行列の最後の要素 A(m,$) - m行の最後の要素 A($,n) - n列の最後の要素 V(m) - ベクトルの第m要素
和 A+B 積 A*B 商 A/B = A*inv(B) 商 A\B = inv(A)*B 行列のべき A^x
成分毎の積 A.*B 成分毎の商 A./B および A.\B 成分毎のべき A.^x
転置複素共役 A' 転置 A.' 複素共役 conj(A)
size(A) - 行列のサイズ[m,n]を与える。 size(A,'r')，size(A,'c')はそれぞれ行列 A の行数，列数を返す． length(A) - 行列の全要素数(すなわちm*n) rank - ランク det - 行列式 trace - トレース cond - 条件数 inv - 逆行列 svd - 特異値分解 sum - 要素の和 cumprod - 累積積 cumsum - 累積和 diag - 対角要素の挿入もしくは抽出 diff - 差分と離散微分 dsearch - 二分探索 find - boolean行列のtrueに対応するindexを与える。 gsort - 順序の並び替え intersect - ２つのベクトルの共通の値のベクトルを返す lex_sort - 辞書式の行列の行の並び替え max - 最大値 maxi - 最大値とそのインデックス min - 最小値 mini - 最小値とそのインデックス nnz - 行列の中の非零の要素の数 gsort - 並び替え
expm - 行列指数関数(exp(A)が要素毎の指数関数であるのに対し、expmは1+A+(1/2)A^2+...を計算する。以下の関数も同様。) acoshm - 行列逆双曲線余弦 acosm - 行列逆余弦 asinhm - 行列逆双曲線正弦 asinm - 行列逆正弦 atanhm - 行列逆正接 atanm - 正方行列逆双曲線正接 coshm - 行列の双曲線余弦 cosm - 行列の余弦 cothm - 行列の双曲線余接 logm - 正方行列の対数 signm - 行列の符号関数 sinhm - 行列の双曲線正弦 sinm - 行列の正弦 sqrtm - 行列の平方根 tanhm - 行列の双曲線正接 tanm - 行列の正接
evals=spec(A)は固有値の配列(列ベクトル)を与え、 [X,B]=spec(A)は固有ベクトルXと固有値を対角成分に持つ行列Bを与える。
poly(m,'x') 配列mを根にもつ多項式を生成する。
例えば
poly([1,2],'x')
は
x^2-3x+2
を生成する。一旦生成された多項式は演算をほどこしても多項式(有理式)のままだから、これは
x=poly(0,'x') p=x^2-3*x+2
としてもよい。
roots - 多項式の根を求める denom - 有理式の分母 numer - 有理式の分子 derivat - 有理式の導関数 horner(p,q) - 有理式pに値qを代入(qは有理式でもよい)
clf - グラフ消去
1 : 黒 2 : 青 3 : 緑 4 : 水色 5 : 赤 6 : マゼンタ 7 : 黄色 8 : 白
-1 : + -2 : x -3 : ○の中に+ -4 : ◆ -5 : ◇ -6 : △ -7 : ▽ -8 : ◇の中に+ -9 : ○ -10 : * -11 : □ -12 : 右向き中空三角 -13 : 左向き中空三角 -14 : ☆
plot2d( x, y [,options] )
x,yは同じ長さの配列。yが行列なら複数のグラフを重ねて描画する。 ただし、scilabのplotは前のグラフを消さないので、無理に行列を 作る必要はなく、単に続けて書けばいい。消すにはclf命令か、ウインドウ のメニューを使う。
○options
style : 各曲線のスタイルをセットする．その値は整数(正あるいは負)の値をもった 実ベクトル． rect : 最小の境界要求をセットする．その値は4つの要素もった実ベクトルである:[xmin,ymin,xmax,ymax] . logflag : 軸のスケール(線形あるいは対数)をセットする． その値はとりうる値の文字列: "nn" , "nl" , "ln" , "ll" . frameflag : 最小の要求値から実際の座標の地域の計算を制御する．その値は0から8 の整数． axesflag : 軸をどのようにするか指定する．その値は0から5の整数． nax : 軸のラベルと目盛りの定義をセットする．その値は4つの整数をもつ実ベクトル [nx,Nx,ny,Ny] leg : 曲線のタイトルをセットする．その値は文字列である．
その他の2次元グラフとして
plot2d2 - ステップ関数として曲線をプロットする plot2d3 - 垂直バーで曲線をプロットする plot2d4 - 矢印で曲線をプロットする fplot2d - ある関数により定義された曲線をプロットする champ - 2次元のベクトル空間 champ1 - 2次元のベクトル空間（色の付いた矢印） fchamp - ある関数により定義された2次元のベクトル空間 contour2d - 2次元平面上への表面の等高線 fcontour2d - 2次元平面上への関数で定義された表面の等高線 grayplot - 色を使った表面の2次元プロット fgrayplot - 色を使った関数で定義された表面の2次元プロット Sgrayplot - 色を使った表面の滑らかな２次元プロット Sfgrayplot - 関数で定義された表面の色を使った滑らかな２次元プロット xgrid - 2次元プロット上にグリッドを加える errbar - 2次元プロット上に縦のエラー・バーを加える histplot - ヒストグラムをプロットする Matplot - 色を使った行列の2次元プロット
がある
例
plot2d(x,y,-9); scatter; plot2d(x,y,2); line; plot2d(x,y,-9,logflag='nl'); log-linear
等高線の例
fname="C030SUM-spectra.dat"; dt=1.72*6; wl0=419.03; dwl=0.28549; m=fscanfMat(fname); x=1239.85./(wl0-dwl*m(:,1)); y=dt*(0:size(m,'c')-2); z=m(:,2:$); clf; contour2d(x,y,z,20);
plot3d - 3次元プロット plot3d1 - 等高線付きの3次元プロット fplot3d - ある関数により定義された3次元プロット fplot3d1 - ある関数により定義された等高線付きの3次元プロット param3d - 曲線をプロットする param3d1 - 曲線を数本プロットする contour - 3次元平面上への等高線 fcontour - 関数で定義された3次元表面上への等高線 hist3d - ヒストグラムの3次元表現 genfac3d - 3次元曲面の面を計算 eval3dp - 3次元パラメトリック曲面の面を計算 geom3d - 3次元プロット後2次元上に3次元からの投影
ワーキングディレクトリに関数ファイルmyfunc.sci等を作成し
getf("myfunc.sci")
あるいは
exec("myfunc.sci")
とすれば,その後使うことができる。myfunc.sciの中身は例えば
function b=myfunc(Xmat,yvec) b=inv(Xmat'*Xmat)*Xmat'*yvec; endfunction
などとなる。オンラインでやるなら
deff("b=myfunc(Xmat,yvec)","b=inv(Xmat'*Xmat)*Xmat'*yvec")
などとする。
戻り値としては変数の組を返すこともでき、その場合の定義は
deff("[a,b]=myfunc(x,y)","a=3*x+2*y, b=4*x-5*y")
のようになる。呼び出し側では、
[aret,bret]=myfunc(x0,y0)
とする。
関数の実行を途中で中断して呼び出し元に戻るには
return;
を使う。あるいは以下で述べるresumeを使っても良い。
グローバルエリアはレベルゼロである。グローバルエリアから呼び出された関数の内部はレベル1、さらに そこから呼び出された関数の内部はレベル2というようにレベルが下がっていく(数が大きい方を下位とする)。 関数内で使用される変数は基本的に全てローカル変数であるが、もし上位レベルに同じ名前の変数があればその値を初期値として引き継ぐ。ただし、その値を関数内で変更しても上位レベルの変数の値には影響を与えない。この機能を使えば、引数を使わなくても関数に値を送ることができる。上位レベルの変数を下位レベルから書き換えるには
[ a1,b1, ... ] = resume( a2, b2, ... )
とする。このとき関数の実行は即座に中断され、上位レベルの変数 a1, b1, ... は a2, b2, ... の内容で書き換えられる。
普通の線形連立方程式 A*x=B は単に x=A\B で解が得られる。 ちなみにA\Bは列ベクトルを与え、B/Aは行ベクトルを与える。
roots(p) は多項式 p について p=0 の根を与える
solve - シンボルの線形システムソルバー
○使い方
x=solve(A,b)
○パラメータ
A : 文字の上三角行列 x,b : 文字の行列
このとき、solveはA*x=bの解を与える。
linsolve – 線形方程式を解く
○使い方
[x0,kerA]=linsolve(A,b [,x0])
○パラメータ
A : n x m の実行列 (場合によりスパース) b : n x 1 ベクトル (Aと同じ行次元) x0 : 実ベクトル kerA : m x k 実行列
linsolveはA*x+b=0の全ての解を計算する。 解は x=x0+kerA*w (wは任意)で与えられる。
(注) 普通の線形連立方程式 A*x=B は単に x=A\B でよい。
fsolve - n変数の非線形関数システムの零点を見つける
○使い方
[x [,v [,info]]]=fsolve(x0,fct [,fjac] [,tol])
○パラメータ
x0 : 列ベクトル(関数引数の初期値)
fct : 目的関数(関数もしくはリストあるいは列)
fjac : fctの導関数(関数もしくはリストあるいは列)
tol : 実数のスカラー。精度誤差 : アルゴリズムが x と解の相対誤差が
大きくとも tol になることを予測すると終了する
( tol=1.d-10 が初期値である )。
○戻り値
x : 実数ベクトル (関数引数の最後の値、およそ零) v : 実数ベクトル (x での関数の値) info : 終了情報 0 : ふさわしくない入力パラメータ 1 : アルゴリズムは x と解の相対誤差が多くとも tol であることを予測する。 2 : fcn に呼び出される数 3 : tol は非常に小さい。近似解ｘのこれ以上の改善は見込まれない。 4 : 繰り返しをしても進歩がない。
fsolveはfct(x)=0となるxをx0の近くに見つける。fjacはfctの導関数を指定する。
(例)
function [y]=myfunc(x) y(1)=sin(x(1))+x(2); y(2)=x(1)-cos(x(2)); endfunction // fsolve([0.7;-0.7],myfunc)
lsq – 線形最小二乗問題
○使い方
X=lsq(A,B)
○パラメータ
A : 実または複素(m x n)行列。mはデータ点の数。nは最適化すべきパラメータの数 B : 実または複素(m x p)行列。pはデータ列の数
○戻り値
X : 実または複素 (n x p) 行列
X=lsq(A,B)は式A*X=Bに対して、ノルム最小となるようなXを 計算する。 例えば f=ax+by+cz+dという関数でデータ点fi(xi,yi,zi)を フィットする場合、Aは
[x1 y1 z1 1]
[x2 y2 z2 1]
[x3 y3 z3 1]
:
[xn yn zn 1]
となり、Bは
[ f1 ] [ f2 ] [ f3 ] : [ fn ]
となる。最適化されたパラメータはXに
[ a ] [ b ] [ c ] [ d ]
の形で与えられる。
複数のデータ列を同時にフィットしたい場合はXもBも複数列とする。
(例)
// // build data // x = (0:0.1:3)'; y = x.*x - x + 3 + (rand(x)-0.5); plot2d(x,y,-2) // //fit // A=[x.*x, x, ones(x)]; B=[y]; X=lsq(A,B); y1=X(1)*x.*x+X(2)*x+X(3); plot2d(x,y1,2)
lsqrsolve - 非線形最小自乗法
○使い方
[p, v] = lsqrsolve(p0, f, m [,fjac])
○パラメータ
f(p,m) : パラメータpについて各データ点での偏差(をベクトルとして)を与える関数 p0 : 初期パラメータの入った列(縦)ベクトル m : データ点の個数
○オプション
fjac : fの導関数
○戻り値
p : 最適化パラメータの列ベクトル v : 列ベクトル。xでの偏差の値。
lsqrsolveはLevenberg-Marquardt法によってsum(f(p,m).^2)を極小にするpを 与える。
(例)
//
// build data
//
a=1;b=0;c=1;
deff('y=FF(x)','y=a*exp(-((x-b)/c)^2)');
X=(-3:.1:3)';
Y=FF(X)+rand(X)-.5;
//
// solve
//
function e=f1(p,m)
a=p(1);b=p(2),c=p(3),
e=Y-a*exp(-((X-b)/c)^2);
endfunction
[p,v]=lsqrsolve([2;1;2],f1,size(X,1));
p
Y2=p(1)*exp(-((X-p(2))/p(3))^2);
plot2d(X,Y,style=-3)
plot2d(X,Y2,style=2)
複数のデータ列を複数の関数でそれぞれ同時にフィットしたい場合には、各データ列の偏差ベクトルをつなげて返すようにf(p,m)を作ればよい。
(例)
//
// Model
//
deff("aa=drude(x,p)","aa=p(1)-p(2)./(x.*x+%i*x./p(3))");
//
// Build data
//
x=(0.5:0.1:10)';
p0=[0, 1.0, 0.5];
y=drude(x,p0);
yr=real(y)+(rand(x)-0.5)/20;
yi=imag(y)+(rand(x)-0.5)/20;
clf;
plot2d(x,yr,-9);
plot2d(x,yi,-11);
//
// Solve
//
function e=err(p,m)
y2=drude(x,p);
e=[yr-real(y2); yi-imag(y2)];
endfunction
p1=[0.1, 0.9, 0.6];
[p,v]=lsqrsolve(p1,err,size(x,1)*2);
//
// Plot
//
y3=drude(x,p);
plot2d(x,real(y3));
plot2d(x,imag(y3));
//
datafit - 離散データへの関数フィッティング
○使い方
[p,err]=datafit( G, Z, p0 )
○パラメータ
G : 関数記述子 (e=G(p,z), e: ne x 1, p: np x 1, z: nz x 1) Z : 行列 [z_1,z_2,...z_n] 。ここで、z_i (nz x 1) は i 番目の測定値である。 p0 : 初期推定 (サイズ np x 1) (オプションについてはマニュアル参照)
○戻り値
p : 列ベクトル。見つけられる最適解。 err : スカラー。最小自乗誤差。
与えられた関数G(p,z)に対して、datafitは一組のベクトルz_iに対して、 G(p,z_i)=0 に近似するパラメータ p の最も良いベクトルを見つける。 アルゴリズムはデフォルトではquasi-Newton法が使われる。
(例)
deff('y=FF(x)','y=a*(x-b)+c*x.*x')
X=[];Y=[];
a=34;b=12;c=14;for x=0:.1:3, Y=[Y,FF(x)+100*(rand()-.5)];X=[X,x];end
Z=[Y;X];
deff('e=G(p,z)','a=p(1),b=p(2),c=p(3),y=z(1),x=z(2),e=y-FF(x)')
[p,err]=datafit(G,Z,[3;5;10])
leastsq - 非線形最小自乗問題を解く
○使い方
[f,xopt]=leastsq(fun [,dfun] ,x0) (その他のオプションは省略)
○パラメータ
fun : 外部関数。すなわち、Scilab 関数もしくは列 ( fun は最小自乗問題を定義している関数である。 dfun : 関数の勾配 x0 : 実数ベクトル (最小化される変数の初期値)。
○戻り値
f : 最適な最小自乗値の値。 xopt : 発見された x の最適値。
制約なしもしくは制約ありのプログラムに対する非線形最適化ルーチン。 fun は関数 f(x) の定義を与える。これは外部関数である。 この外部関数は与えられた x に対して y(j)=f(x)のような y を返さなければならない。 アルゴリズムはデフォルトではquasi-Newton法が使われる。
(例)
a=rand(3,2);b=[1;1;1];x0=[1;-1];
deff('f=fun(x,a,b)','f=a*x-b');
deff('g=dfun(x,a,b)','g=a');
[f,xopt]=leastsq(fun,x0) //Simplest call
xopt-a\b //compare with linear algebra solution
[f,xopt]=leastsq(fun,dfun,x0) //specify gradient
[f,xopt]=leastsq(list(fun,[1 2;3 4],[1;2]),x0)
deff('f=fun(x,a,b)','f=exp(a*x)-b');
deff('g=dfun(x,a,b)','g=a.*(exp(a*x)*ones(1,size(a,2)))');
[f,xopt]=leastsq(list(fun,[1 2;3 4],[1;2]),x0)
制約なしの関数最小化。Nelder–Mead法を使ってパラメータを最適化し、コスト関数を最小化する。
○使い方
xopt = fminsearch ( fun , x0 )
○パラメータ
fun : コスト関数 x0 : 実数ベクトル(パラメータの初期値)
○戻り値
xopt : 発見されたパラメータの最適値
interp - 3次元のスプライン評価関数 interp2d - 2次元のスプライン評価関数 interp3d - 3次元のスプライン評価関数 interpln - 線形補間 lsq_splin - 加重最小自乗法の３次元スプラインフィッティング bsplin3val - 関数の任意の微分係数の3次元スプライン cshep2d - 2次元 3次 shepard (scattered)補間 eval_cshep2d - 3次元の shepard 補間計算 linear_interpn - n次元の線形補間 splin - 3次元スプライン補間 splin2d - スプラインの格子付き２次元補完 splin3d - スプラインの格子付き３次元補完
integrate - 求積法による積分 intsplin - スプライン補間によるデータの積分 inttrap - 台形補間によるデータの積分 delip - 楕円積分 int2d - 求積法もしくは立体求積法による２次元の定積分 int3d - 求積法もしくは立体求積法による３次元の定積分 intc - コーシー積分 intg - 定積分 intl - コーシー積分
ode - 常微分方程式ソルバ ode_discrete - 常微分方程式ソルバー、離散時間シミュレーション ode_root - 根を求める常微分方程式 odedc - 離散/ 連続時間 ode ソルバー odeoptions - ode ソルバーのオプションを決める
微分方程式 dy/dx=f(x,y)を解く。
使い方
ode(y0,x0,x,f);
y0 : yの初期値。連立微分方程式の場合は縦ベクトル。 x0 : xの初期値。スカラー。 x : 求めたい結果のx座標。ベクトル。 f : 関数。連立微分方程式の場合は縦ベクトル関数。
odeはx0から適当なステップで数値的に微分方程式を解き、 xで与えられたx座標に対するyの値を返す。xのステップと 計算時のステップは異なる。
(例)1階常微分方程式
deff("ydot=f(x,y)","ydot=-y");
x0=0;
y0=1;
x=x0:0.1:10;
y=ode(x0,y0,x,f);
plot2d(x,y);
(例2)1階連立常微分方程式
deff("ydot=f(x,y)","ydot=[-3*y(1)-4*y(2);-2*y(1)-5*y(2)]");
x0=0;
y0=[0;1];
x=0:0.1:10;
y=ode(y0,x0,x,f);
plot2d(x,y');
1階の連立常微分方程式に変形して解けばいい。
