Microsoft PowerPoint - it2_octave_text1.ppt

Transcription

1 数値計算ソフトウェア 背景 機械技術者は様々な力学現象を扱う. この現象を表わす解が要求された場合, 先ずは解析解 ( 厳密解 を求める. しかし, 一般には困難である. 解析解が得られない場合には, 工学的に許容される数値近似解を求めることになる. 目標様々な数値計算を行ない, 汎用数値計算ソフトウェアを活用できるようにする. 1 数値計算数値近似解 行列演算 非線型方程式 常微分方程式 統計解析 制御系設計グラフ表示可視化 gnuplotなどの外部ソフトウェアを用いて表示プログラムバッチ処理 if,for,whileなどの制御構造 代表的な汎用数値計算ソフトウェア と Maxima の特徴的な相違点 処理方式 解 値の型 グラフ 備考 数値計算 数値近似解 実数型 数値データ 複雑なシステム大規模なシステム Maxima 数式処理 解析解 ( 厳密解 整数型実数型 関数 性質の良い問題 商用 Matlab 非商用 Scilab 豊富な関数群 ビジュアル化による直感的理解 高度な演算処理 応用数値計算に関する専用ソフトの例 気象予測, 衝突解析, 応力解析, 流体解析, 伝熱解析 ショートカットの作成 ( 4 プロパティの変更 ( (1 octave のショートカット を右クリック 5 の起動 6 (1 のショートカット を選択 ( 送る ( 作業フォルダ の指定を変更 (4 デスクトップ ( ショートカット %HMEDRIVE%%HMEPATH%octave フォルダ z: octave が作成された ( 初回のみ表示される ( を右クリック moodle に文字列が記載されているので, 複写 貼り付けするとよい コマンドの入力場所 (1 スタートメニュー すべてのプログラム edsys アプリケーション GNU ( K を選択 プロンプト 紙面の都合で以降は >> と記載する の終了 算術演算子 終了 quit octaveの終了 exit octaveの終了 ヘルプ help コマンド等の一覧 help command 指定 commandの解説表示 7 demo プログラムの実行 demo ファイルを実行確認しなさい. 1 moodle の指定場所から, 圧縮ファイル octave_sample.zip を Download し, 展開する. 8 ファイル octave_demo.m example.m を, フォルダ z: octave に保存する. で octave_demo と入力し, 実行する. 基本的な処理 9 ( 注意 Window の強制終了ではエラーが発生する. quit あるいは exit により終了すること. >> octave_demo スクリプトファイル octave_demo.m 処理を記述したテキストファイル 操作方法 演算子 関数

2 入力方法の例 結果 : 入力出力 操作方法 (peration (+5*7/-1 ans を計算に利用できる >> ans * ans 79 >> (+5*7/-1 7 一時記憶変数 ; は不要 入力した計算式結果 大文字と小文字を区別 10 算術演算子 (Arithmetic operator 算術演算子例結果 + 和 x+y 差 x-y 5- * 積 x*y 5* 15 / 商 x/y 5/ ^ 冪 x^y 5^ 15 ** 冪 x**y 5** 15 は, 整数, 小数などの数値を全て倍精度浮動小数点の形式 ( 近似値 で記憶, 計算する 注意 : 打切り誤差, 丸め誤差に注意計算方法により結果が異なることがある 11 数学定数 (Constants 定数 値, 備考 pi 円周率 e ネイピア数 ( 自然対数の底 i,j,i,j 虚数単位 0+1i eps 計算機の精度 e-16 realmax 記憶可能な最大実数 e+08 realmin 記憶可能な最小実数 e-08 Maxima とは異なり, pi, e に % を付けない. 1 数学関数 (Mathematical Functions 1 三角関数 (Trigonometric Functions 14 計算例結果 abs(x 絶対値 abs(- sqrt(x 平方根 sqrt(^+4^ 5 exp(x 指数関数 exp( log(x 自然対数 log(e^ log10(x 常用対数 log10(100 ceil(x 切り上げ ceil(.5 sign(x 符号 sign(- -1 引数 戻り値 sin(x 正弦 x[rad] -1~1 cos(x 余弦 x[rad] -1~1 tan(x 正接 x[rad] inf~-inf asin(x 逆正弦 -1~1 -pi/~pi/ acos(x 逆余弦 -1~1 0~pi atan(x 逆正接 ( 象限 inf~-inf -pi/~pi/ atan(y,x 逆正接 (4 象限 inf~-inf x atan(x 1 y x atan(y,x -pi~pi 複素数 (Complex Number, 複素関数 計算例結果 a+b*i,a+bi 複素数 +4*i +4i real(z 実部 real(+4*i imag(z 虚部 imag(+4*i 4 abs(z 大きさ abs(+4i 5 arg(z 偏角 arg(1+i conj(z 共役 conj(+4*i -4i 大きさ sqrt(a^+b^ 偏角 atan(b,a b Im 偏角 a a+b*i Re 共役 a-b*i 15 定義方法 変数の定義 (= 変数名 = 値あるいは計算式 ( 注 変数は値を代入しなければ使えない. Maxima との特徴的な相違点. >> z error: `z' undefined... 例 変数の削除 (clear >> x=5 >> y= y = >> x+y 変数 z には値が記憶されていないため, エラー発生 指定した変数を削除 clear 変数全ての変数を削除 clear 16 数値リスト ベクトル型 行列型 数値リストの演算 数値リスト 縦ベクトル型 行列型 横ベクトル 行列 [x1,x,x,..,xn] [x1;x;x;..;xn] [a11,a1,..,a1n; a1,a,..,an;...; am1,am,..,amn] >> x=[5,,7,4] x = >> A=[1,,;4,5,6] A = 区切り, 要素 ; 行 空白 改行 縦ベクトル 要素 行 >> x=[4;;;6] x = 4 6

3 数値リストの生成 ( 増加量指定 備考初期値 : 増加量 : 範囲 範囲内で打ち切り 第 項に減少量を指定しても良い 処理例 初期値 : 範囲増加量は 1 >> 1:0.4: >> 0: 数値リストの生成 ( 分割数指定 分割形式備考等間隔 linspace(a,b,n 区間 [a,b] を n 個に等分割. n の省略時は n=100. 対数 logspace(a,b,n 区間 [10^a,10^b] について [a,b] を n 個に等分割. ( 注 両端を含めて n 個 >> linspace(0,pi, 数値リストの生成 ( 行列型 要素の値正方行列 m n 行列例 (m=,n= 全要素が 0 zeros(n zeros(m,n [0,0,0; 0,0,0] 対角要素が1 eye(n eye(m,n [1,0,0; 0,1,0] 全要素が 1 ones(n ones(m,n [1,1,1; 1,1,1] 対角行列 diag([x1,x,..,xn] [x1, 0,.., 0; 0,x,.., 0;...; 0, 0,..,xn] 1 演算子 数式 + 和 A+B - 差 A-B * 積 A*B * スカラーとの積 p*a ^ 冪 A^p ' 複素共役転置 A'.' 転置 A.' M zm1 z1n M z 行列の演算 転置 共役転置 A,B: 行列 p: スカラー 複素数 共役複素数 M z1n M z1n z= x+ yi z = x yi zm1 M z zm1 M z 行列の演算関数 ( 長方形 計算結果転置 transpose(a [1,4;,5;,6] 行列のサイズ size(a [,] 行のサイズ rows(a 列のサイズ columns(a 階数 rank(a A=[1,,; 4,5,6] 行列式 逆行列 対角要素の和 固有値分解 行列の演算関数 ( 正方行列 A=[1,;,4]; [v,lambda]=eig(a A= det(a inv(a trace(a [v,lambda]=eig(a [[ v1, vn ],diag[ λ1,, λn ]] 4 固有値と固有ベクトル v= lambda= A A a 11+ a+ + a nn Avi = λivi ( i= 1,,, n 線型連立方程式次の線型連立方程式の解を求めなさい. x+ y= x+ 4y= 7 行列を用いて, 係数と未知数を分離する. x = 4 y 7 逆行列を用いて解を求める x = = y 4 7 / 1 = / Ax=b x=a b ベクトルの追加 拡大係数行列 リストの積, 追加 処理 式 結果 ベクトルの積 u*v' u'*v [ 4, 5, 6; 8,10,1; 1,15,18] [u,w] [1,,,4,5,6] [u;w] [A,b] [1,,; 4,5,6] [1,,5;,4,6] u=[1,,] v=[4,5,6] A=[1,;,4] b=[5;6] 6 演算子.* 積./ 商.^ 冪 f( 関数 要素毎の演算 C=A./B C=A.^B C=f(A 計算例 >> x=[,4,5,6], y=[8,6,4,] x = y = >> x.* y 数式計算結果 C=A.*B cij=aij*bij cij=aij/bij cij=aij^bij cij=f(aij A,B は同サイズの数値リスト 7 A=[1,;,4];b=[;7]; x=inv(a*b x = >> x./ y

4 A=[,4;5,6]; B=[8,6;4,]; 要素毎の関数演算 >> sqrt(a >> sin(b 要素毎の関数演算 A= 要素毎の積 8 B= 4 6 >> sqrt(a.* sin(b ( 参考 通常の行列の積 >> sqrt(a * sin(b グラフ 次元曲線 (plot 次元曲線 (plot 次元網目 (mesh 9 次元グラフ (plot 0 次元平面内に曲線を描画する plot([x1,..,xn],[y1,..,yn] x=-10:0.1:10; plot(x,sin(x; y= sinx ( 0 x 次元グラフ ( 重ねて描画 plot(x_list,[y1_list;..;ym_list] x=-10:0.1:10; plot(x,[sin(x;cos(x;sin(*x+cos(*x]; sinx y= cosx sin x+ cosx ( 0 x 1 次元グラフ (plot 次元空間内に曲線を描画する plot([x1,..,xn],[y1,..,yn],[z1,..,zn] t=linspace(0,10,00; plot(t.*sin(pi*t,t.*cos(pi*t,t; ( x, y, z = ( t sinπt, t cosπt, t (0 t 要素毎の積に注意 次元グラフ (mesh 網目状に描画する mesh( [x1,..,xn],[y1,..,ym], [z11,..,z1n;..;zm1,..,zmn] x=y=linspace(0,*pi,50; sx=sin(x; cy=cos(y; z=cy' * sx; mesh(x,y,z; z= sinx cosy (0 x π,0 y π [y1,.., ym] 座標の [x1,..,xn] [z11,..,z1n;...; zm1,..,zmn] z ij = f ( xj, yi 作業フォルダ 4 例 path 探索フォルダの一覧 addpath(path 探索フォルダへpathを追加 rmpath(path 探索フォルダからpathを削除 genpath(path pathのサブフォルダを抽出 dir 作業フォルダのファイル一覧 cd 作業フォルダを ユーザーフォルダへ変更 cd dir 作業フォルダをdirへ変更 拡張子 スクリプトファイル 5 注意点 m 保存場所現在の作業フォルダ, あるいは path で一覧されるフォルダ 記載内容対話形式で入力する手順と同じ内容. ただし, 最初にダミー命令を記載する. 実行方法 上でファイル名を入力する. スクリプトファイルの例 example.m 1; x = 0:0.1:10; y = sin(x; plot(x,y ダミー命令 対話形式で入力する手順と同じ disp(arg 画面表示 input('msg' キーボード入力 デフォルトの表示桁数 表示桁数の変更 内部の演算桁数 表示設定を戻す 入出力関数 使用方法 disp(x disp('msg' a=input('msg' 表示桁数 (format 値 x の表示 文字列の表示 値の入力 >> pi pi =.1416 >> format long >> pi pi = >> format 有効数字は 15 桁程度 6

5 >> x=5 >> y= y = >> x+y 命令の区切り 区切り結果の画面表示 改行 次の行に分ける場合 する, 1 行に続ける場合 する ; 1 行に続ける場合しない 改行の場合, の場合 ; の場合 >> x=5,y=,x+y y = 各処理の結果が画面表示される >> x=5; >> y=; >> x+y; >> x+5;y=;x+y; 画面表示が抑制される 7 管理者対応 S John W. Eaton inux, Windows, MacS,... 記述言語 C++, Standard Template ibrary (ST ライセンス GNU GP 1988 年頃化学反応装置に関する学部用教科書と共に使用するソフトウェアを意図して,James B. Rawlings, Jhon W. Eaton により作成された 199 年 より柔軟な数値計算ツールとして本格的に開発 1994 年 Ver1.0がリリースされる 1997 年 現在 GNU GP に移行し, 自由に使えるようになる. 現実的な問題を解くために, 多様な分野の講義で利用されている 8 参考文献 9 (1 John W. Easton: GNU Manual ( Quick Reference ( 北本 : を用いた数値計算入門, ピアソン エデュケーション (4 吉田和信 : による動的システムシミュレーション入門, ctave.htm (5 松田七美男 : の概要, s.dendai.ac.jp/~matuda/tex/lecture.html (6 北野正雄 : 関数電卓とのすすめ, kyoto-u.ac.jp/jp/common/course4/pdf/octave.pd f