×

[PR]この広告は3ヶ月以上更新がないため表示されています。
ホームページを更新後24時間以内に表示されなくなります。

■Excel:ソルバーの使い方・・・方程式,行列,固有値問題などを解く
 ここでは,Excelに付いている「ソルバー」を使って,方程式,行列計算,固有値問題などを解くことを考える.・・・マクロ(プログラム)は使わない.
 以下に述べる解説は,それぞれの項目のおよその意味が分かることを前提とする.もし分からない項目があれば,そこに立ち返るか,その項目は省略するとよい.

 たとえば,

(1) 2次方程式の解き方は分かるが,「 x2 + 109.67x - 1325.6 = 0 を筆算で解くのは大変なのでExcelにやらせよう」という場合

(2) 「行列の固有値問題とは何のことかだいたい分かるが,実数解のうち大きい方から2,3個求めたい」という場合

 のように,数学の問題を数値的に解くことを考える.
※ 「ソルバー」を使うためのExcelの設定

○ メニュー→「ツール」で「ソルバー」が表示されていればOK.

▲ メニュー→「ツール」で「ソルバー」が表示されていなければ,「ツール」→「アドイン」で「ソルバーアドイン」にチェックを付ける.

(これもできないときは,OfficeのCD-ROMを使ってインストールから始める.)
■1 ソルバーの使い方に慣れよう
■例1 2次関数の最小値の求め方
 2次関数 y = ax2 + bx + c はa>0のとき右下図のようなグラフとなり,最小値がある.適当にa,b,cの値を入力して,「ソルバー」を使って最小値を求めたい.

 y = 1x2 + 2x + 3 の最小値を求めたいとき,

(1) 右のようにx,a,b,cの値を入力しておく.(この場合xは初期値)

(2) セルF2には2次関数の計算式を書き込んでおく:
=B2*A2^2+C2*A2+D2

(3) メニューから「ツール」→「ソルバー」→
目的セルをF2とする
目標値で「最小値」を選択
変化させるセルをA2とする
実行ボタン」をクリック
解を記入する」でOK
  A B C D E F
1 x= a= b= c=   ax^2+bx+cの値:
ソルバーで目標
とするセル
2 0 1 2 3   3
3 ↑ソルバーで変化
させるセル
         

(参考)
 y = x2 + 2x + 3 =(x + 1)2 + 2 は,
x=-1 のとき最小値 y=2 をとる.
「ソルバー」を何度も使うとき,前のパラメータ(ユーザが指定できる値)がソルバーの画面に残っているので,「ツール」→「ソルバー」→「リセット」からスタートするとよい.

■例2 2次方程式の解の求め方

 1x2 + 2x - 3 = 0 の解を求めたいとき,

(1) 右のようにx,a,b,cの値を入力しておく.(この場合xは初期値)

(2) セルF2には2次関数の計算式を書き込んでおく:
=B2*A2^2+C2*A2+D2

(3) メニューから「ツール」→「ソルバー」→
目的セルをF2とする
目標値で「」 0 を選択
変化させるセルをA2とする
実行ボタン」をクリック
解を記入する」でOK

※ 初期値(A2の値)を十分大きな値(10など)にして再度実行すると大きい方の解x=1が得られる.
  A B C D E F
1 x= a= b= c=   ax^2+bx+cの値:
ソルバーで目標
とするセル
2 -10 1 2 -3   77.0
3 ↑ソルバーで変化
させるセル
         

(参考)
 2次方程式の「実数解」は,ないことも,1個のことも,2個のこともある.

x2 + 2x - 3 = 0の解は
(x + 3)(x - 1) = 0より
x = -3,1 の2つある.
初期値x = -10からスタートすればx=-3が得られ,
初期値x = 10からスタートすればx=1が得られる.(上図赤の矢印参照)
 Excelのソルバーでは,変化させるセルの値を初期値から少しずつ変化させて,目的セルの値が指定された値になれば止まる.「解なし」のときは,無限ループとはならず「仮の解が見つかりません」と表示される.
■例3 条件式を減らす工夫

(1) 0に等しいという条件
x = 0,y = 0 ⇔ x2 + y2 = 0 
x = 0,y = 0,z = 0 ⇔ x2 + y2 + z2 = 0
・・・
のように2つ以上の条件式は1つの式に直すことができる.
(2) 各々別の値に等しいという条件
x = a,y = b ⇔ (x-a)2 + (y-b)2 = 0 
x = a,y = b,z = c ⇔ (x-a)2 + (y-b)2 + (z-c)2 = 0
・・・
のように2つ以上の条件式は1つの式に直すことができる.
(2)の例として,連立方程式 ax + by = c, dx + ey = f の解をソルバーを利用して求めるには

(1) 右のようにA3:E4にa,b,c,d,e,f,x,yの値を入力しておく.(この場合x,yは初期値)

(2) セルG3:G4にはax+by-c,dx+ey-fの計算式を書き込んでおきH3:H4にはそれらの2乗を,H5に縦の和を書き込んでおく.

(3) メニューから「ツール」→「ソルバー」→
目的セルをH5とする
目標値で「」 0 を選択
変化させるセルをE3:E4とする
実行ボタン」をクリック
解を記入する」でOK
  A B C D E F G H
1 a= b= c=   x=   ax+by-c 2乗
2 d= e= f=   y=   dx+ey-f   
3 1 -2 5   2.999999903   0.0000 0.0000
4 4 5 7   -1.000000073   0.0000 0.0000
5             0.0000

※ この問題の解は,数学的にはx = 3,y = -1であるが,上の例のように x = 2.9999・・・,y = -1.0000・・・のようになることがある.小数の数値計算では,このように厳密な値でなく所定の誤差の範囲内にあるものが解となる.
 誤差の範囲は「オプション」で設定を変えれば調整できるが,実測値中心の次のような計算では,精度は重要であるが表記は気にならない.

1.4142 -1.7321 = 2.2361
2.4495 + 2.6417 = 2.828
→  = 1.3542
   = -0.1852
■例4 目的セルを複数個指定したいとき
「ソルバー」では,
目的セル」は1個しか指定できないが,
変化させるセル」は「1続きの範囲」「飛び飛びのセルをカンマで区切って」いずれも複数個指定できる
制約条件」は複数個指定できる

 複数個の条件を指定する例として,ソルバーを使って逆行列を求めることを考えてみる・・・なお,逆行列は,関数MINVERSE()でも求められるので,結果が合っているかどうか点検できる.(※逆行列の意味についてはこのページ参照)

 行列Aに対して,AB=E(単位行列)となる行列BをAの逆行列という.
 右図のように元の行列Aの成分をA2:B3に書く.
 逆行列Bの成分を書き込む場所をD2:E3とする(ただし,初期値はでたらめに書いておく)
 行列の積ABの成分をH2:I3に書く.(*)
 (点検用の行列はD6:E7などに書いておく.(**)

(1) 右のように行列Aの成分を入力し,次にその逆行列Bの成分はでたらめに入れておく.(これが初期値)

(2) セルH2:I3には行列ABの積の計算式を書き込んでおく.

(3)
 BがAの逆行列であるためには,
H2=1,I2=0,H3=0,I3=1
という形で4個の条件式が必要
 しかし,目的セルは1つしか指定できない.これをどうするか?
==> 1つは目的セルとして指定,残り3個は制約条件にすれば複数個指定可能.(各々の差の2乗の和を作れば1つにまとめることも可能)

 メニューから「ツール」→「ソルバー」→
目的セルをF3とする
目標値で「」 0 を選択
変化させるセルをA2とする
制約条件で残り3個のセルの値を指定する
実行ボタン」をクリック
解を記入する」でOK
  A B C D E F G H I
1 A=元の行列を適当に指定   B=ソルバーで変化させるセル   AB=「目標とするセル」「制約条件」に指定する
2 1 2   -5 2     1.000 0.000
3 3 5   3 -1     0.000 1.000
4                  
5       点検用:MINVERSE()
関数による解
         
6       -5 2        
7       3 -1        
(*)
行列の積の求め方
H2に=MMULT(A2:B3,D2:E3)と書き込むと(1,1)成分だけが書き込まれるので,次にH2:I3を反転表示させてから,上端の数式バーにマウスを当て,Ctrl+Shift+Enterとする.(ややむずかしい)
(**)
逆行列を関数で求める方法
D6に=MINVERSE(A2:B3)と書き込むと(1,1)成分だけが書き込まれるので,次にD6:E7を反転表示させてから,上端の数式バーにマウスを当て,Ctrl+Shift+Enterとする.(ややむずかしい)

※ 作者のExcelは,ソルバーで制約条件を「追加」→「OK」とするとエラー表示(入力した参照が正しくない・・・)となるが,さらに「キャンセル」→「実行」で動く.
■2 ソルバーで統計の問題を解いてみよう
■ソルバーでは「最小値」を求めることが」できるので,「近似直線・近似平面」を実際に「最小2乗法」によって直接求めることができる.
 回帰直線の方程式を求める問題 ( correl.xls )において,右のように与えられたx,yの値の組(40組)について最小2乗法を用いて近似直線(回帰直線)の方程式を求めるには,
(1) 求める直線の方程式を y = ax+ b とおき,a,bの値を書き込むセルを作り,適当に初期値を入れておく(右図ではA44:B44)
(2) 各観測値xに対してyの予測値 ax+b を求める.
D2に=a$44*a2+b$44 とし,これを下までコピー・貼り付けする.
(3) 求めた予測値とyの観測値との誤差を計算する.
E2に=B2-D2 と記入し,これを下までコピー・貼り付けする.
(4) 誤差の2乗を求め、次にその和を求める.
F2に=E2^2 と記入し,これを下までコピー・貼り付けする.
E44に=SUM(E2:E41)とする.

 メニューから「ツール」→「ソルバー」→
目的セルをF44とする
目標値で「最小値」 を選択
変化させるセルをA44:B44とする
実行ボタン」をクリック
解を記入する」でOK
  A B C D E F
1 x y   ax+b y-(ax+b) その2乗
2 13.031 21.969   21.786 0.183 0.033
3 12.968 25.291   21.679 3.611 13.041
4 14.610 19.506   24.466 -4.960 24.603
・・・ ・・・ ・・・   ・・・ ・・・ ・・・
42            
43 a= b=      
44 1 1       3711.36 
■ソルバーでは「最大値」を求めることができるので「主成分分析」の考え方で分散が最大となる係数の組合せ(第一主成分)を直接求めることができる.
 主成分分析では,複数の変数に座標変換を行って,分散が最大となるような新変数(主成分)を作る.
 資料間の差異が一番はっきり分かる変数(主成分)を作る = 分散が最大となる新変数を求める =
X = ax + by とおいてXの分散が最大となる係数a,bを定める.

 これは図のように回転させて新変数X の分散が最大となるように変形することに対応している.

 右のように2つの観測値をもつ資料が50個あって,x,y単独では資料間の差異がはっきりしないとき,
X = ax + by となる新変数Xを導入して,その分散が最大となるようにする.
(1) 求める変数をX = ax + by とおき,a,bの初期値を右のようにB54:C54にでたらめな数字で書き込む.
(2) それらの2乗と2乗の和をB56:D56に計算式で書く.
 B56には =B54^2,C56には =C54^2
 D56には =SUM(B56:C56)
(3) 各資料についてax+byの計算式をE列に書く.
 E2に =B$56*B2+C$56*C2
この式をコピーしてE3からE51までに貼り付け.
(4) E列の分散の計算式をE56に書く.
 =VAR(E2:E51)

 メニューから「ツール」→「ソルバー」→
目的セルをE56とする
目標値で「最大値」を選択
変化させるセルをB54:C54とする
制約条件で$D$56 = 1 とする.
実行ボタン」をクリック
解を記入する」でOK
        ===> (第一主成分が求まる.)
(高校数学の復習)
X = ax + by
Y = cx + dy
の関係式によって古い点(x,y)に新しい点(X,Y)を対応させるものを1次変換という.
 今の高校の数学では,(x,y)→(X,Y)を点の移動として教えるが,左図のように点がそのままで,座標軸を動かすと考えても同じになる.その場合,「各点を原点の周りに30°回転させること」は「座標軸を原点の周りに-30°回転させること」と同じになる. [1次変換の解説][回転移動の解説

 2次元で原点の周りの回転のとき次が成り立つ. 
X = ax + by (a2 + b 2 = 1)
 3次元で原点の周りの回転のとき次が成り立つ. 
X = ax + by + cz (a2 + b 2 + c2 = 1)
 n次元で原点の周りの回転のとき次が成り立つ. 
X1 = a1x1 + a2x2 + a3x3 + ・・・ + anxn
  (a12 + a2 2 + a32 + ・・・ + an2= 1)

  A B C D E
1   x y   ax+by
2 資料1 -5.1 14.04   13.369
3 資料2 16.9 20.04   21.832
4 ・・・ ・・・ ・・・   ・・・
51 資料50 5.9 -9.96   -9.224
53   a= b=    
54   1 1    
55    a^2= b^2= 分散
56    1 1 2 5682
別添ファイル (solver1.xls)のSheet1
■ソルバーでは目的セルの「」を指定できるので,固有値問題を解くことができる.

○1 固有値・固有ベクトルを定義に従って直接求めることもできる・・・ただし,固有ベクトルは大きさ自由
Ax=λx → A(kx)=λ(kx)
だから大きさ1の条件を付ける.
Ax=λx (ただし,|x| = 1)

(1) 行列A(右の表1では2×2行列)の成分を入力しておく.(A2:B3)
(2) 固有ベクトルxの初期値を適当に入力しておく.
(C2:C3)
(3) 固有値の初期値を適当に入力しておく.(D2)
(4) Ax(A6:A7)およびλx(B6:B7)の計算式を書き込む.これらの差Ax-λxの計算式も書き込む.(C6:C7)
(5) 固有ベクトルの大きさ(の2乗)を計算しておく.
 F2に =C2^2,F3 に=C3^2,F5に =SUM(F2:F3)

 メニューから「ツール」→「ソルバー」→
目的セルをC6とする
目標値で「」 0 を選択
変化させるセルをC2:C3 カンマを打ってさらにD2
制約条件で$F$5 = 1,C7 = 0 とする.
実行ボタン」をクリック
解を記入する」でOK


○2 特性方程式(固有方程式)を解いて固有値を求めることができる:行列Aの固有方程式 | A - λE | = 0 の解を求める.
(1) 行列A(右の表2では2×2行列)の成分を入力しておく.(A2:B3)
(2) 固有値の初期値を適当に入力しておく.(D2)
(4) λE(A5:B6)の計算式を書き込む.差A-λEの計算式も書き込む.(A8:B9)
(5) | A-λE | を計算しておく.
 D9 に=MDETERM(A8:B9)

 メニューから「ツール」→「ソルバー」→
目的セルをD9とする
目標値で「」 0 を選択
変化させるセルをD2とする.
実行ボタン」をクリック
解を記入する」でOK
表1
  A B C D E F
1 A=   x= λ=    xの2乗
2 3 6 1.000 10    1.000
3 2 -1 1.000      1.000
4          
5 Ax= λx= AX-λx=     2.000
6 9.000 10.000 -1.000       
7 1.000 10.000 -10.000        

※ 行列  の固有値はλ=-3,5
2つとも見つけるには
初期値を10としてスタートすれば,最大の固有値5が得られる.
次に,D7の制約条件4.9以下として追加,再度実行すると-3も得られる.







表2
  A B C D
1 A=     λ=
2 -2 5   5.000
3 4 -3    
4 λE=      
5 5.000 0.000    
6 0.000 5.000    
7 A-λE=      
8 -7.000 5.000   |A-λE|=
9 4.000 -8.000   36.000

※ 行列  の固有値はλ=-7,2
2つとも見つけるには
上と同様に求めた値よりも少し小さい値以下を制約条件にすれば順次求まる.
■例と答■
 
(1) 次の方程式の実数解を1つ見つけよ.
11x5 - 10x4 + 9x3 - 8x2 + 7x - 6 = 0
(答)
x = 0.884
(2) 次の方程式の実数解を1つ見つけよ.
sinx + sin2x + sin3x = cosx cos2x cos3x
(参考)
y=sinx+sin2x+sin3x-cosx*cos2x*cos3xのグラフは周期2πの周期関数となり
区間0≦x<2πでは,0.1461, 1.5708(=π/2), 2.2013, 2.8478, 4.3811, 4.7124(=3π/2)が方程式を満たす.
(3) 次の行列の固有値を4個求めよ.
0 2 -2 1
-3 5 1 -2
-2 1 0 2
1 -2 -3 5
(答)
λ = 6,3,1,0 の4個あるが
初期値10からスタートしてλ=6
次に,制約条件λのセルの値に <= 5.9 とするとλ=3
次に,制約条件λのセルの値に <= 2.9 とするとλ=1
次に,制約条件λのセルの値に <= 0.9 とするとλ=0
が得られる.
(4) 次の資料(別添ファイル solver1.xls のSheet2)について,
X = ax + by + cz + du +ev
(a2 + b2 + c2 + d2 + e2= 1)
となる新変数Xを導入して,その分散が最大となるように定数 a,b,c,d,eを求めよ.
 さらに,ベクトル (a,b,c,d,e)に垂直なベクトル(p,q,r,s,t)について
Y = px + qy + rz +su + tv
(p2 + q2 + r2 + s2 + t2 = 1)
となる新変数Yを導入して,その分散が最大となるように定数 p,q,r,s,tを求めよ.
(データは http://www.jma.go.jp/jma/index.html 気象庁>過去の気象データ検索>京都市の月ごとの気温,湿度等の2004〜2006年を標準化したもの)
(答)
変数 x y z u v
  降水量
(mm)
気温
(℃)
湿度
(%)
風速
(m/s)
日照時間
(h)
X 係数 a b c d e
  -0.175 0.196 -0.574 0.515 0.580
Y 係数 p q r s t
  0.639 0.645 0.248 0.330 -0.072

 このXYで散布図を作ると次のようになる.
 月の区別は温度でできると簡単に考えていたが,そうではなく,第1主成分のXは風速,日照時間-湿度に特徴があり,第2主成分のYは温度,降水量に特徴がある.

○=== メニューに戻る

◇このページの内容について,考え方の間違い,計算間違い,著作権上の
問題点などお気づきの点がございましたら までご連絡ください.