【第5回】

<内容>

1.        数学関数

·       sincosなどの基本的な数学関数が用意されている.これらの関数について学ぶ.

·       数学関数を用いた簡単な関数のグラフ表示の方法を学ぶ.

 

2.        多項式に関する関数

·       多項式の表現方法や,零点,特性多項式の求め方を学ぶ.

·       与えられたデータから最小二乗近似した多項式曲線の求め方を学ぶ.

·       多項式の最大公約数の求め方,ならびに部分分数展開の方法を学ぶ.

 

 

<演習>

1.数学関数

代表的な数学関数が2.6(教科書p.38)に示されている.例えば,y = sin(t)のグラフを表示させる場合は,以下のようにする.

  t = 0:0.1:10;        % 0.1刻みで,010を要素とするベクトルtが生成される.

  y = sin(t);       % sin(t)の値を要素とするベクトルyが生成される.

  plot(t,y)         %横軸をt,縦軸をyとするグラフを描く.

  grid(’on’)        % グラフに格子を描く.デフォルトでは’off’

 

(注)plotコマンドを最初に実行した際,グラフが描画されるまでに時間が掛かる.2回目以降は,すぐに描画される.

 

 

グラフウィンドウをクリアするには「clf」(clear figure)を,グラフウィンドウを閉じるには,「close」を入力する.

 

  テキストでは,plotgridの順番がすべて逆になっている.また,Octave 4.0では,closeplotの代わりにcloseを使用する.

 

  Octave 4.0では,plotを実行しても凡例(line 1line 2等)がグラフの右上に表示されない.凡例を表示させるには,plotの後にlegendを用いる.

 

(例)線が1つの場合:legend('line 1')
線が2つの場合:legend('line 1','line 2')

line 1line 2がグラフに表示される.これらの文字は自由に変更可能.

 

 

練習5-1

 例題2.10(教科書p.38),例題2.11(教科書p.39)を行え.また,例題2.11nを大きくしていくと波形がどのようになっていくかを確認せよ.

 

注)y = exp(-0.5*t).* cos(5*t)の掛け算「.*」に注意すること.exp(-0.5*t)cos(5*t)ベクトル(1×1001である(下表参照).

掛け算を単に「*」とするとベクトル(行列)同士の掛け算を意味し,次元が合わないので計算できない.演算子の前に「.」を置くと,それは要素ごとの演算を意味し,対応する要素ごとに掛け算が実行される.

 

t

exp(-0.5*t)

cos(5*t)

0.00

1.00000

1.00000

0.01

0.99501

0.99875

0.02

0.99005

0.99500

10.00

0.00674

0.96497

 

 

問題5-1

 問題1.(教科書p.40)を行え.また,nを大きくしていくと波形がどのようになっていくかを確認せよ.(ヒント:例題2.11

 

注)ベクトル(行列)とスカラーとの演算では,常に要素ごとの演算となる点に注意すること.Matlabでも同様である.

(例)

> t = [0  1  2];

> 1 – t

ans =

   1  0  -1

> 2 * t

ans =

   0  2  4

> A = [ 1  2 ; 3  4 ];

> A – 1

ans =

   0  1

   2  3

 

 

2.多項式に関する関数

 多項式に関する代表的な関数が2.7(教科書p.42)に示されている.各関数の説明が見たい場合は,「help 関数名」を入力すればよい.例えば,関数rootsの場合は

  help  roots

となる.ただし,説明文はすべて英語である.

 

 

21.多項式の表現,零点,特性多項式

 多項式は,その係数を降べき順(次数が高い→低い)に並べたベクトルとして表される.例えば,

  p(x) = -x4 + 20x2 - 20x + 5

を表すには

  p = [ -1  0  20  -20  5 ];

とすればよい.この多項式の零点p(x) = 0の根)は関数「roots」を用いて次のように計算できる.

  r = roots(p)

実行結果は以下のようになり,ベクトルrに4つの零点が入力される.

r =

   -4.92606

    3.89858

    0.57356

    0.45393

逆に,零点から多項式の係数を求めるには

  p1 = poly(r)

とする.実行結果は

  p1 =

      1.0000e+000  3.4972e-015  -2.0000e+001  2.0000e+001  -5.0000e+000

となる.ベクトルp1に多項式の係数が降べき順に入力される.ただし,最高次数の係数を「1」としたものになる点に注意せよ.ここでは,p1-pとなっている.

 

 

練習5-2

 上記のことを確認してみよ.また,roots(p)roots(p1)の値が同じであることも確認せよ.

 

 

 関数「poly」(polynomial:多項式)は正方行列A特性方程式を求めるときにも使用できる.例えば,

  図面1.jpg

の場合

  A = [ 0  1  0 ; -1  0  0 ; 0  -1  2];

  p2 = poly(A)

とすればよい.行列Aの特性多項式の定義式は

  図面1.jpg

である.ここで,| | は行列式を意味し,I は単位行列を示す.行列A固有値eigenvalue)は多項式p2(x)の零点であるので,

  eig(A)

  roots(p2)

は同じ値となる.すなわち,eig(A) = roots(poly(A))である.

 

 

練習5-3

 上記のeig(A)roots(poly(A))が同じ値となることを確認してみよ.

 

 

 多項式p(x)x = 1における値は,関数「polyval」(polynomial value)を用いて

  y = polyval(p,1)

で計算できる.

 

 

練習5-4

 例題2.12(教科書p.43)を行え.

 

 

―――――――――――――――――――――――――――――――――――――――――――――――――――――――

 

 

22.多項式の曲線のあてはめ

 2次元データから最小二乗近似を用いてn次多項式を求める関数が「polyfit」(polynomial fitting)である.

 

 

練習5-5

次のデータ

  x = [ 0  1  2  3  4  5  6  7  8  9  10 ]

  y = [ 0  3  4  10  18  27  41  59  79  93  101 ]

n次多項式曲線で近似するプログラム「ex5.m」(教科書p.44)を作成し,n = 1, 2, 3とした場合の結果を比較してみよ.

 

【補足】

同一グラフ上にn本の線を描く場合は

  plot(x1,y1,x2,y2, ,xn,yn)

のようにすればよい.また,各x, yの組のあとに '+' 'o' のオプションを書くと線を引かずに,+ o がプロットされる.

 

 

23.多項式の最大公約数,部分分数展開

 多項式 f(x)g(x) ユークリッドの互除法

   f = gq1 + r1                 %  f / g = q1 余り r1

g = r1q2 + r2                %  g / r1 = q1 余り r2

r1 = r2q3 + r3               %  r1 / r2 = q3 余り r3

     …

rn-1 = rnqn+1                %  rn-1 / rn = qn+1 余り 0

 を適用すると,( f(x)g(x) ) の最大公約数GCDgreatest common divisor)がrnとして得られる.

 

 

練習5-6

 例題2.13(教科書p.45)を行え.

 

 <解説>

·    多項式の掛け算の関数が「conv」(convolution:畳み込み),多項式の割り算の関数が「deconv」(deconvolution:反畳み込み)である.

·    [q, r] = deconv(f, g) で,f/gを行い,商がqに,剰余がrに代入される.

·    fgは多項式の係数ベクトルであるので,[2  3  2  1][1 3 4 2]のように入力する.

·    Matlabには,高次の余分な0要素を除いた多項式を求める関数「polyreduce」が用意されていない.関数polyreduceを使用しない場合を以下に記しておく.これならMatlabでも実行可能である.

f=input('Enter f ');

g=input('Enter g ');

while(1)       % 無限ループ:break文が実行されるまで永遠にループする

    [q,r]=deconv(f,g);

    if norm(r)<1e-10

        break;

    end

    while(abs(r(1))<1e-10)           % if文でなくwhile文を用いる点に注目せよ

        r(1)=[]; % 先頭行から0の要素が全て削除される

    end

    f=g;

    g=r;

end

gcd=g/g(1)

 

 

 有理関数の部分分数展開を行う関数「residue」(留数)は,多項式 P(s)Q(s) に対して,具体的には下記の各係数rpkeを出力する.

  図面1.jpg

 

 

練習5-7

 例題2.14(教科書p.46)を行え.特に,residue(P,Q)の結果から最終的な解

  図面1.jpg

が得られることを理解せよ.ただし,k=[](0x0)は,k00列の空行列であることを意味する.

  実行結果の表示順がテキストと異なるが,実質的には全く同じものである.

 

 

問題5-2

(1)      問題1.(教科書p.47)を行え.(ヒント:根を求めれば,因数分解できたことになる.)

(2)      問題2.(教科書p.47)を行え.(ヒント:例題2.13

(3)      問題3.(教科書p.48)を行え.(ヒント:例題2.14

 

【答え】

(1)     

(2)      図面1.jpg

(3)     

 

 

【レポート課題(5)(授業終了後にアップロード)

  RePORT