アルゴリズム演習2 MATLABプログラミング

Presentation on theme: "アルゴリズム演習2 MATLABプログラミング"— Presentation transcript:

1 アルゴリズム演習2 MATLABプログラミング
岩崎　慶

2 演習内容 MATLABプログラミングの習得 スクリプトMファイル ファンクションMファイル 非線形連立方程式の解法

3 Mファイル Ｍファイルには2種類存在する スクリプトMファイル ファンクションMファイル 入出力変数をもたない、コマンドを集めたファイル
入出力変数をもつ、関数を定義したファイル 関数名とファイル名は同一である必要がある サーチパス内にあれば別の関数からも呼び出すことが可能 「ファイル」→「パスの設定」で作成した関数を保存したディレクトリをパスに含める 今回は簡単のため、自分の適当なディレクトリ内にファンクションＭファイルを作成し、そのディレクトリ内で使用する

4 Mファイルの例 スクリプトMファイル ファンクションMファイル 出力変数 入力変数 関数名(ファイル名と同じ)
x = 0:pi/100:2*pi; y = sin(x); plot(x,y); A = [4 0 -1; 0 5 0; ]; b = [1; 2; 0]; x = A \ b function f = func( x, y, z ) f = x + y + z; func.m function f = func( x, y, z ) 出力変数 入力変数 関数名(ファイル名と同じ) test.m >> func(1,2,3) ans = 6

5 ファンクションMファイル 関数ごとにプログラムを作成するのは無駄 𝑓と𝑔が違うだけで構造は同じ 関数𝑓とその導関数𝑓1 初期値𝑥0を
f x^2-2; f1 2*x; x = 2; %f(x)が0になるまで計算 while(abs(f(x)) > 0.001) x = x – f(x) / f1(x); end g x^3+ 4; g1 3*x^2; x = 2; %g(x)が0になるまで計算 while(abs(g(x)) > 0.001) x = x – g(x) / g1(x); end 関数ごとにプログラムを作成するのは無駄 𝑓と𝑔が違うだけで構造は同じ 関数𝑓とその導関数𝑓1 初期値𝑥0を 与えたら解を返す関数を作れないか？

6 ファンクションMファイル 出力変数 引数(入力変数) 使用法 f = @(x) x^2-2; f1 = @(x) 2*x;
while(abs(f(x)) > 0.001) x = x – f(x) / f1(x); end function x = newton(f, f1, x0) x = x0; %f(x)が0になるまで計算 while(abs(f(x)) > 0.001) x = x – f(x) / f1(x); end newton.m 関数名とファイル名は同じ 使用法 f x^2 – 2; f1 2*x; x = newton(f, f1, 2); 引数に渡す関数を変えれば様々な関数で ニュートン法による解が求められる

7 無名関数のように関数を渡すことができない
ファンクションMファイル function y = func_g(x) y = 2*x^2 + 3 * log(x) – 1; function y = func_g1(x) y = 4*x + 3 / x; function x = newton(f, f1, x0) x = x0; while(abs(f(x)) > ) x = x – f(x) / f1(x); end newton( func_g, func_g1, 2 ) エラー 無名関数のように関数を渡すことができない

8 関数ハンドル 関数の計算に必要な関数の全てのデータをもつMATLABのデータタイプ 関数を他の関数の引数に渡すことができる @関数名
関数ハンドルの評価 g g1 g( 1 ); feval(g, 1);

9 ニュートン法による非線形連立方程式の解法
例　2変数( 𝑥 1 , 𝑥 2 )連立方程式の場合 初期解候補 𝒙 (0) と近傍の値 𝒙 (0) +𝜹𝒙を以下のように定める 2変数関数の1次のテーラー展開を 𝒙 (0) の近傍点で行う 𝑓 1 𝑥 1 , 𝑥 2 =0, 𝑓 2 𝑥 1 , 𝑥 2 =0 𝒙 (0) = 𝑥 1 (0) , 𝑥 2 (0) , 　𝛿 𝒙 (0) =(𝛿 𝑥 ,𝛿 𝑥 2 (0) ) 𝑓 1 𝒙 0 +𝛿 𝒙 0 = 𝑓 1 𝒙 0 + 𝜕 𝑓 1 𝒙 (0) 𝜕 𝑥 1 𝛿 𝑥 1 (0) + 𝜕 𝑓 1 𝒙 (0) 𝜕 𝑥 2 𝛿 𝑥 2 (0) 𝑓 2 𝒙 0 +𝛿 𝒙 0 = 𝑓 2 𝒙 0 + 𝜕 𝑓 2 𝒙 (0) 𝜕 𝑥 1 𝛿 𝑥 1 (0) + 𝜕 𝑓 2 𝒙 (0) 𝜕 𝑥 2 𝛿 𝑥 2 (0)

10 非線形連立方程式 𝑓 1 ( 𝒙 (0) +𝛿 𝒙 (0) ) 𝑓 2 ( 𝒙 (0) +𝛿 𝒙 (0) ) = 𝑓 1 ( 𝒙 (0) ) 𝑓 2 ( 𝒙 (0) ) + 𝜕 𝑓 1 𝒙 (0) 𝜕 𝑥 1 𝜕 𝑓 1 𝒙 (0) 𝜕 𝑥 2 𝜕 𝑓 2 𝒙 (0) 𝜕 𝑥 1 𝜕 𝑓 2 𝒙 (0) 𝜕 𝑥 𝛿 𝑥 1 (0) 𝛿 𝑥 2 (0) ヤコビ行列 0 0 = 𝑓 1 ( 𝒙 (0) ) 𝑓 2 ( 𝒙 (0) ) + 𝜕 𝑓 1 𝒙 (0) 𝜕 𝑥 1 𝜕 𝑓 1 𝒙 (0) 𝜕 𝑥 2 𝜕 𝑓 2 𝒙 (0) 𝜕 𝑥 1 𝜕 𝑓 2 𝒙 (0) 𝜕 𝑥 𝛿 𝑥 1 (0) 𝛿 𝑥 2 (0)

11 非線形連立方程式 𝜕 𝑓 1 𝒙 (0) 𝜕 𝑥 1 𝜕 𝑓 1 𝒙 (0) 𝜕 𝑥 2 𝜕 𝑓 2 𝒙 (0) 𝜕 𝑥 1 𝜕 𝑓 2 𝒙 (0) 𝜕 𝑥 𝛿 𝑥 1 (0) 𝛿 𝑥 2 (0) = − 𝑓 1 ( 𝒙 0 ) − 𝑓 2 ( 𝒙 0 ) この式により𝛿 𝒙 (0) を求め, 𝒙 (1) = 𝒙 (0) +𝜹 𝒙 (0) と解に近づける 𝛿 𝒙 (𝑘) <𝜀もしくは 𝒇 𝒙 𝑘 =( 𝑓 1 𝒙 𝑘 𝑓 2 𝒙 𝑘 ) <𝜀となったら反復終了

12 非線形連立方程式 非線形連立方程式を解くプログラム newton_m.m
% f: 関数ハンドル df:ヤコビ行列の関数ハンドル x0:初期値　 function x = newton_m( f, df, x0 ) x = x0; %初期値代入 count = 0; %ループのカウント数 dx = - df( x ) \ f( x ); %𝛿𝑥を計算 while( norm( dx ) > ) x = x + dx; dx = - df( x ) \ f( x ); count = count + 1; if( count > 100 ) break; %無限ループに陥らないため強制終了 end

13 非線形連立方程式 𝑥 1 2 + 𝑥 2 2 − 3 =0 3 𝑥 1 2 +7 𝑥 2 2 −8=0 を解く
𝑥 𝑥 2 2 − 3 =0 3 𝑥 𝑥 2 2 −8=0 　を解く gfunction.m function f = gfunction(x) f = [ x( 1 )^2 +x( 2 )^2-sqrt( 3 ); 3*x( 1 )^2+7*x( 2 )^2-8 ]; d_gfunction.m function f = d_gfunction(x) f = [ 2 * x( 1 ) 2 * x( 2 ); 6 * x( 1 ) 14 * x( 2 ) ]; g dg newton_m( g, dg, [3;4])

14 disp関数 文字列または配列の表示 disp(‘this is a test’); %文字列の表示 x = 1; %変数(1行1列の行列)
disp(x); A = [1 2; 3 4]; disp(A); %行列の表示 disp(‘ ‘);%改行の表示(クォーテーションの間はスペース)

15 sprintf関数 定数や変数を文字列に変換する関数 変換指定子 例 変数・定数 sはd文字列 実行結果
s = sprintf(‘one degree is %f \n’, pi / 180); disp(s); sはd文字列 実行結果 one degree is

16 fprintf関数 定数や文字列を画面に表示する関数 例 結果
fprintf( ‘one degree is %f\n’, pi / 180 ); >> fprintf( 'one degree is %f\n', pi / 180 ); one degree is 結果

17 変換指定子 意味 %c 単一キャラクタ %d int型の10進数表現（符号付） %e
指数表現　%Eは3.1415E+00のように大文字Eの表示 %f 固定小数表示　%7.5fは全桁数7桁、小数点5桁 %g %eまたは%fのどちらか短い表現 %s 文字列キャラクタ表現 %o 8進数表現（符号なし） %u 10進数表現（符号なし） %x 16進数表現

18 ヤコビ法 𝐴=𝐷+(𝐿+𝑈)と分解 𝐷 : 対角行列 （対角成分以外0） 𝐿 : (狭義)下三角行列 𝑈 : (狭義)上三角行列
𝐷 : 対角行列　（対角成分以外0） 𝐿 : (狭義)下三角行列 𝑈 : (狭義)上三角行列 𝑈 𝐷 𝐿 𝐷= 𝑎 11 ⋯ 0 ⋮ ⋱ ⋮ 0 ⋯ 𝑎 𝑛𝑛 𝐷 −1 = 1/𝑎 11 ⋯ 0 ⋮ ⋱ ⋮ 0 ⋯ 1/𝑎 𝑛𝑛 𝐿= 0 ⋯ 0 ⋮ ⋱ ⋮ 𝑎 𝑛1 ⋯ 0 𝑈= 0 ⋯ 𝑎 1𝑛 ⋮ ⋱ ⋮ 0 ⋯ 0

19 ヤコビ法 𝒙 (𝑘+1) = 𝐷 −1 {𝒃− 𝐿+𝑈 𝒙 (𝑘) } 𝒙 (𝑘+1) = 𝐷 −1 𝒃− 𝐿+𝑈+𝐷 𝒙 𝑘 +𝐷 𝒙 𝑘
𝒙 (𝑘+1) = 𝐷 −1 {𝒃− 𝐿+𝑈 𝒙 (𝑘) } 𝒙 (𝑘+1) = 𝐷 −1 𝒃− 𝐿+𝑈+𝐷 𝒙 𝑘 +𝐷 𝒙 𝑘 𝒙 (𝑘+1) = 𝐷 −1 𝒃−𝐴 𝒙 𝑘 + 𝐷 −1 𝐷 𝒙 𝑘 𝐴 𝒙 (𝑘+1) = 𝐷 −1 𝒃−𝐴 𝒙 𝑘 + 𝒙 𝑘

20 ヤコビ法のヒント 𝒙 (𝑘+1) = 𝐷 −1 (𝒃−𝐴 𝒙 𝑘 )+ 𝒙 𝑘 𝒓=𝒃−𝐴 𝒙 𝑘
𝒙 (𝑘+1) = 𝐷 −1 (𝒃−𝐴 𝒙 𝑘 )+ 𝒙 𝑘 𝒓=𝒃−𝐴 𝒙 𝑘 𝒙 (𝑘+1) は 𝒙 𝑘 に 𝐷 −1 𝒓を足して更新 𝒓 が閾値以下になったら反復終了 %A : 係数行列, b : 右辺のベクトル, threshold : 閾値 %D : Aの対角成分からなるベクトル matrixsize = size( A, 1 ); x = zeros( matrixsize, 1 ); maxiteration = 100; epsilon = zeros( 1, maxiteration ); for k = 1 : maxiteration % rをA, x, bから計算 epsilon( k ) = norm( r ); %norm( r )が閾値thresholdより小さければ終了 %xをD, rから計算し，x^(k+1)に更新 end