Maple と MATLABを使った振動解析

Maple は, 数式処理, 数値計算, グラフィック, プログラムの機能を統合したソフトウェアです。教育機関での教育・研究から企業・研究機関での開発・設計まで幅広い場面で利用されています。

数式処理ソフトウェアを使用するメリットとしては,

等の点が挙げられます。これらの理由から, 純粋に数学的な問題を追及する用途以外にも, 理論式の把握や確認, 数値計算ソフトウェアに渡す前の下処理など, 研究開発の基本ツールとして活用することができます。

また, MATLAB リンクの機能を使用することで, Mapleで作成した数値行列を MATLABのカーネルを使用して計算できるようになりました。
この機能を利用すると, モデル式の作成からMATLABでの数値解析までの作業を, 1つのインターフェース上で行えるようになります。
ここでは, MATLABリンク機能を利用して多自由度振動の解析を行う例をご紹介します。

>    restart;

1. Mapleを使ったモデル式の作成

[Maple OLE 2.0 Object]
図1

n個のおもりをばねでつなぎ, 両端を固定したときの自由振動を考えます。(図1) おもりの質量はすべて m とします。

例えば n=3のとき, Mapleではこのシステムを次のように記述できます。ここでは xの2階微分を簡単のため ddxと書きます。

>    [k1*x1-k2*(x2-x1)=-m*ddx1,
k2*(x2-x1)+k3*(x2-x3)=-m*ddx2,
k3*(x3-x2)+k4*x3=-m*ddx3];

[k1*x1-k2*(x2-x1) = -m*ddx1, k2*(x2-x1)+k3*(x2-x3) = -m*ddx2, k3*(x3-x2)+k4*x3 = -m*ddx3]

Mapleの特徴は, 変数の入った数式をそのまま扱えることです。ここでは, 先ほどの式を, collect コマンドを使って x1, x2, x3 についてまとめます。

>    collect(%,[x1,x2,x3]);

[(k1+k2)*x1-k2*x2 = -m*ddx1, -k2*x1+(k2+k3)*x2-k3*x3 = -m*ddx2, -k3*x2+(k3+k4)*x3 = -m*ddx3]

次に, おもりの数がn個の場合の運動方程式は次のように表すことができます。
 
      
(k1+k2)*x1-k2*x2 = m*ddx1   ----   n = 1  
      
-(k.n*x.(n-1))+(k.n+k.(n+1))*x.n-(k.(n+1)*x.(n+1)) = m*ddx.n   
      
-(k.N*x.(N-1))+(k.N+k.(N-1))*x.N = m*ddx.N     ----   n = N

式1
ここでは, Mapleのプログラミング機能を使い, おもりがN個のときの連立方程式を生成するプロシージャmakesys を作ってみましょう。

>    makesys:=proc(N,de)
 local j;
 de:=array(1..N);
 de[1]:=(k1+k2)*x1-k2*x2 = -m*ddx1;
 for j from 2 to N-1 do
   de[j]:=-k||j*x||(j-1)+(k||j+k||(j+1))*x||j-k||(j+1)*x||(j+1) = -m*ddx||j;
 od:
   de[N]:=-k||N*x||(N-1)+(k||N+k||(N+1))*x||N = -m*ddx||N;
   print(de);
 RETURN();
end:

例えば, n=5では, 次の式が生成され, 変数deに代入されます。

>    makesys(5,de);

vector([(k1+k2)*x1-k2*x2 = -m*ddx1, -k2*x1+(k2+k3)*x2-k3*x3 = -m*ddx2, -k3*x2+(k3+k4)*x3-k4*x4 = -m*ddx3, -k4*x3+(k4+k5)*x4-k5*x5 = -m*ddx4, -k5*x4+(k5+k6)*x5 = -m*ddx5])
vector([(k1+k2)*x1-k2*x2 = -m*ddx1, -k2*x1+(k2+k3)*x2-k3*x3 = -m*ddx2, -k3*x2+(k3+k4)*x3-k4*x4 = -m*ddx3, -k4*x3+(k4+k5)*x4-k5*x5 = -m*ddx4, -k5*x4+(k5+k6)*x5 = -m*ddx5])

>   

2. 剛性行列の生成

この式は線形なので, 行列とベクトルを使い

[K] x= -[M] ddx  - A

と表わせます。[K] は剛性行列, [M] は質量行列と呼ばれます。また x, ddxは各質点の変位と加速度を表すベクトルになります。

Mapleの線形代数パッケージには, 線形の連立方程式から係数行列を取り出すコマンドgemnatrixが用意されています。これを使って 行列 [K]と[M]を取り出してみましょう。

まず, 先ほど makesysで作った方程式deより, [K] を取り出します。

>    with(linalg):

Warning, the protected names norm and trace have been redefined and unprotected

>    K:=genmatrix(convert(de,list),[x1,x2,x3,x4,x5]);

K := matrix([[k1+k2, -k2, 0, 0, 0], [-k2, k2+k3, -k3, 0, 0], [0, -k3, k3+k4, -k4, 0], [0, 0, -k4, k4+k5, -k5], [0, 0, 0, -k5, k5+k6]])

同様に右辺の [ mddx1, mddx2, ... mddx5]から, M を取り出してみましょう。

>    M:=genmatrix([m*ddx1, m*ddx2, m*ddx3, m*ddx4, m*ddx5],[ddx1,ddx2,ddx3,ddx4,ddx5]);

M := matrix([[m, 0, 0, 0, 0], [0, m, 0, 0, 0], [0, 0, m, 0, 0], [0, 0, 0, m, 0], [0, 0, 0, 0, m]])

最後に, このシステムの特性をMATLABで求めるため, [K]と[M]に適当な値を代入し, 数値行列を作ります。

ここでは k1.. k6は 0.5, m=1.0とします。

>    K1:=subs({k1=0.5,k2=0.5,k3=0.5,k4=0.5,k5=0.5,k6=0.5},
     evalm(K));

K1 := matrix([[1.0, -.5, 0, 0, 0], [-.5, 1.0, -.5, 0, 0], [0, -.5, 1.0, -.5, 0], [0, 0, -.5, 1.0, -.5], [0, 0, 0, -.5, 1.0]])

>    M:=subs(m=1.0,evalm(M));

M := matrix([[1.0, 0, 0, 0, 0], [0, 1.0, 0, 0, 0], [0, 0, 1.0, 0, 0], [0, 0, 0, 1.0, 0], [0, 0, 0, 0, 1.0]])

>   

3. MATLABリンク機能

MATLABリンクを使うと, Mapleの数値行列を, 同じマシン上のMATLABで計算することができます。直接, 数値計算ルーチンにアクセスするため,特に大きな行列でもMaple 上で高速に解を得ることができます。(精度は MATLABに依存します。)

リンク用のコマンドは, Mapleの「Matlab関数パッケージ」に収録されています。このパッケージを使うには, あらかじめwith(Matlab): コマンドを実行するか, Matlab[コマンド名]のようにパッケージ名を明記してコマンド入力を行います。

コマンドは大きく分けて次の2種類があります。

@ ハイレベルコマンド
MATLABのeig, det, fft, qr, lu, chol, ode45 などのコマンドは, Maple のインターフェース上から直接計算させることができます。
 
[Maple OLE 2.0 Object]

図2

A ローレベルコマンド
MapleとMATLABのワークスペース間で数値行列の
やり取りを行ったり, 指定したMATLABコマンドを実
行させることができます。


[Maple OLE 2.0 Object]

                 図3

4. モード行列の計算

ここで, 式2において x = a*exp(I*t*w)と置き換えます。ここで a は各おもりの振幅を表します。また I は虚数単位を表します。

すると 式2は, [K] a = w2* [M] aとなり, 左から M^2*I をかけ, 変形すると

([M]-I*[K] + w2*I)* a = 0  - 式3

となります。

各おもりが振幅a で振動するには, w が det(M-I*K-w^2*I)=0 - 式4
を満している必要があります。

式3, 4より, w, a  はそれぞれ, [M]-I * [K] の固有値(固有周波数)と固有ベクトル(モードベクトル)に相当していることが分かります。これは, 自由振動が固有周波数においてのみおこることを意味しています。

では,さっそくMATLABリンク機能を利用してこのおもりの様子を調べてみましょう。

>    with(Matlab):  # MATLABリンクの開始

Warning, the previous bindings of the names det and transpose have been removed and they now have an assigned value

>    setvar("K1",K1);

>    setvar("M",M);

>    evalM("[P,W]=eig(K1,M)");

>    P:=getvar("P");

P := Matrix(%id = 37052980)
P := Matrix(%id = 37052980)
P := Matrix(%id = 37052980)
P := Matrix(%id = 37052980)
P := Matrix(%id = 37052980)
P := Matrix(%id = 37052980)

>    W:=getvar("W"):

これでモード行列 [P]と固有周波数 [W]が計算できました。
また, Mapleのコマンドを使っても, 固有値と固有ベクトルを求めることができます。解は,

[固有値, 対応する固有ベクトルの数, {固有ベクトル}]

の形式で表示されます。

>    eigenvects(inverse(M)&*K1);

[1.866025404, 1, {vector([.2886751347, -.4999999997, .5773502689, -.4999999999, .2886751348])}], [.1339745968, 1, {vector([-.2886751343, -.4999999999, -.5773502690, -.5000000003, -.2886751345])}], [1.5...
[1.866025404, 1, {vector([.2886751347, -.4999999997, .5773502689, -.4999999999, .2886751348])}], [.1339745968, 1, {vector([-.2886751343, -.4999999999, -.5773502690, -.5000000003, -.2886751345])}], [1.5...
[1.866025404, 1, {vector([.2886751347, -.4999999997, .5773502689, -.4999999999, .2886751348])}], [.1339745968, 1, {vector([-.2886751343, -.4999999999, -.5773502690, -.5000000003, -.2886751345])}], [1.5...
[1.866025404, 1, {vector([.2886751347, -.4999999997, .5773502689, -.4999999999, .2886751348])}], [.1339745968, 1, {vector([-.2886751343, -.4999999999, -.5773502690, -.5000000003, -.2886751345])}], [1.5...
[1.866025404, 1, {vector([.2886751347, -.4999999997, .5773502689, -.4999999999, .2886751348])}], [.1339745968, 1, {vector([-.2886751343, -.4999999999, -.5773502690, -.5000000003, -.2886751345])}], [1.5...

MapleとMATLABリンクの計算結果を比較してみて下さい。

次に, MATLABリンクの結果から, 式2の中のベクトル a を取り出します。

ここでは, モードベクトルを v1.. v5, それに対応する固有周波数をw1..w5として取り出します。

>    for j from 1 to 5 do
   v||j:=col(P,j): W||j:=W[j,j]:
end do;

v1 := vector([-.2886751346, -.5000000000, -.5773502692, -.5000000000, -.2886751346])

W1 := .133974596215561264

v2 := vector([-.5000000000, -.5000000000, .3243377990e-15, .5000000000, .5000000000])

W2 := .499999999999999944

v3 := vector([.5773502692, -.9464193699e-16, -.5773502692, -.1513948001e-15, .5773502692])

W3 := 1.

v4 := vector([-.5000000000, .5000000000, -.2974833629e-16, -.5000000000, .5000000000])

W4 := 1.50000000000000000

v5 := vector([.2886751346, -.5000000000, .5773502692, -.5000000000, .2886751346])

W5 := 1.86602540378443882

これらのベクトルをMapleのlistplotコマンドでグラ
フにすると, 各おもりの変位を見ることができます。
(図4はv1のグラフ)

>    with(plots):

Warning, the name changecoords has been redefined

>    listplot(v1,view=[0..6,-1..1],title=cat("w=",convert(W1,string)),
style=point,symbol=circle);

[Maple Plot]

図4

また, 次のコマンドを入力すると, おもりの動きをアニメーションで見ることができます。

>    vp1:=evalm(v2*sin(2*Pi/10*theta));

vp1 := vector([-.5000000000*sin(1/5*Pi*theta), -.5000000000*sin(1/5*Pi*theta), .3243377990e-15*sin(1/5*Pi*theta), .5000000000*sin(1/5*Pi*theta), .5000000000*sin(1/5*Pi*theta)])

>    v1[1];

-.2886751346

>    seq(plot([[10+vp1[1],0],[20+vp1[2],0],[30+vp1[3],0],[40+vp1[4],0],[50+vp1[5],0]],style=point,symbol=circle,axes=boxed),theta=1..10):

>    display(%,insequence=true);

[Maple Plot]

5.  n=15 の場合

おもりの数が増えても, 同じ方法でモードベクトルと固有周波数を得ることができます。
n=15のときのモードベクトルの一部を下に示します。

>    makesys(15,de15);

vector([(k1+k2)*x1-k2*x2 = -m*ddx1, -k2*x1+(k2+k3)*x2-k3*x3 = -m*ddx2, -k3*x2+(k3+k4)*x3-k4*x4 = -m*ddx3, -k4*x3+(k4+k5)*x4-k5*x5 = -m*ddx4, -k5*x4+(k5+k6)*x5-k6*x6 = -m*ddx5, -k6*x5+(k6+k7)*x6-k7*x7 =...
vector([(k1+k2)*x1-k2*x2 = -m*ddx1, -k2*x1+(k2+k3)*x2-k3*x3 = -m*ddx2, -k3*x2+(k3+k4)*x3-k4*x4 = -m*ddx3, -k4*x3+(k4+k5)*x4-k5*x5 = -m*ddx4, -k5*x4+(k5+k6)*x5-k6*x6 = -m*ddx5, -k6*x5+(k6+k7)*x6-k7*x7 =...
vector([(k1+k2)*x1-k2*x2 = -m*ddx1, -k2*x1+(k2+k3)*x2-k3*x3 = -m*ddx2, -k3*x2+(k3+k4)*x3-k4*x4 = -m*ddx3, -k4*x3+(k4+k5)*x4-k5*x5 = -m*ddx4, -k5*x4+(k5+k6)*x5-k6*x6 = -m*ddx5, -k6*x5+(k6+k7)*x6-k7*x7 =...
vector([(k1+k2)*x1-k2*x2 = -m*ddx1, -k2*x1+(k2+k3)*x2-k3*x3 = -m*ddx2, -k3*x2+(k3+k4)*x3-k4*x4 = -m*ddx3, -k4*x3+(k4+k5)*x4-k5*x5 = -m*ddx4, -k5*x4+(k5+k6)*x5-k6*x6 = -m*ddx5, -k6*x5+(k6+k7)*x6-k7*x7 =...
vector([(k1+k2)*x1-k2*x2 = -m*ddx1, -k2*x1+(k2+k3)*x2-k3*x3 = -m*ddx2, -k3*x2+(k3+k4)*x3-k4*x4 = -m*ddx3, -k4*x3+(k4+k5)*x4-k5*x5 = -m*ddx4, -k5*x4+(k5+k6)*x5-k6*x6 = -m*ddx5, -k6*x5+(k6+k7)*x6-k7*x7 =...
vector([(k1+k2)*x1-k2*x2 = -m*ddx1, -k2*x1+(k2+k3)*x2-k3*x3 = -m*ddx2, -k3*x2+(k3+k4)*x3-k4*x4 = -m*ddx3, -k4*x3+(k4+k5)*x4-k5*x5 = -m*ddx4, -k5*x4+(k5+k6)*x5-k6*x6 = -m*ddx5, -k6*x5+(k6+k7)*x6-k7*x7 =...

>    kitei15:=[seq(x||j,j=1..15)];

kitei15 := [x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x13, x14, x15]

>    K:=genmatrix(convert(de15,list), kitei15);

K := matrix([[k1+k2, -k2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [-k2, k2+k3, -k3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, -k3, k3+k4, -k4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, -k4, k4+k5, -k5, 0...

>    K15:=subs({k1=0.5,k2=0.5,k3=0.5,k4=0.5,k5=0.5,k6=0.5,k7=0.5,k8=0.5,k9=0.5,k10=0.5,k11=0.5,k12=0.5,k13=0.5,k14=0.5,k15=0.5,k16=0.5},
     evalm(K));

K15 := matrix([[1.0, -.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [-.5, 1.0, -.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, -.5, 1.0, -.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, -.5, 1.0, -.5, 0, 0, 0...

>    M15:=band([1.0],15);

M15 := matrix([[1.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 1.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 1.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 1.0, 0, 0, 0, 0, 0, 0, 0, 0, ...

>    setvar("K15",K15);

>    setvar("M15",M15);

>    evalM("[P,W]=eig(K15,M15)");

>    P:=getvar("P");

P := Matrix(%id = 88891892)

>    W:=getvar("W");

W := Matrix(%id = 89749260)

>    for j from 1 to 15 do
  v||j:=col(P,j): W||j:=W[j,j]:
end do:

>    listplot(v1,view=[0..15,-1..1],title=cat("w=",convert(W1,string)),
style=point,symbol=circle);

[Maple Plot]

>    listplot(v2,view=[0..15,-1..1],title=cat("w=",convert(W2,string)),
style=point,symbol=circle);

[Maple Plot]

6. 応用 - 建物の振動

応用例として, 理想化されたN階建ての建物の例を考えます。(図5)
建物はどのような揺れを起こすでしょうか。

[Maple OLE 2.0 Object]

床に各階の質量が集中し, 壁がばねの働きをすると考えると, 先程のおもりの例と同様に方程式を立てることができます。ただし, 最上階が固定されていないため, 式が若干異なり, 次のようなプロシージャで作成できます。

>    makesys2:=proc(N,de)
  local j;
  de:=array(1..N);
  de[1]:=2*k*x1-k*x2 = -m*ddx1;
  for j from 2 to N-1 do
    de[j]:=-k*x||(j-1)+ 2*k*x||j-k*x||(j+1) = -m*ddx||j;
  end do:
  de[N]:=-k*x||(N-1)+k*x||N = -m*ddx||N;
  print(de);
  RETURN();
end proc:

では, さっそく連立方程式を作成し, この建物の剛性行列を求めてみましょう。ここではk=1.25e8, m=5000として計算します。また, 出力が多いため, コマンドの最後を:(コロン)にし, 結果を表示しないようにします。

>    makesys2(22,de22);

vector([2*k*x1-k*x2 = -m*ddx1, -k*x1+2*k*x2-k*x3 = -m*ddx2, -k*x2+2*k*x3-k*x4 = -m*ddx3, -k*x3+2*k*x4-k*x5 = -m*ddx4, -k*x4+2*k*x5-k*x6 = -m*ddx5, -k*x5+2*k*x6-k*x7 = -m*ddx6, -k*x6+2*k*x7-k*x8 = -m*dd...
vector([2*k*x1-k*x2 = -m*ddx1, -k*x1+2*k*x2-k*x3 = -m*ddx2, -k*x2+2*k*x3-k*x4 = -m*ddx3, -k*x3+2*k*x4-k*x5 = -m*ddx4, -k*x4+2*k*x5-k*x6 = -m*ddx5, -k*x5+2*k*x6-k*x7 = -m*ddx6, -k*x6+2*k*x7-k*x8 = -m*dd...
vector([2*k*x1-k*x2 = -m*ddx1, -k*x1+2*k*x2-k*x3 = -m*ddx2, -k*x2+2*k*x3-k*x4 = -m*ddx3, -k*x3+2*k*x4-k*x5 = -m*ddx4, -k*x4+2*k*x5-k*x6 = -m*ddx5, -k*x5+2*k*x6-k*x7 = -m*ddx6, -k*x6+2*k*x7-k*x8 = -m*dd...
vector([2*k*x1-k*x2 = -m*ddx1, -k*x1+2*k*x2-k*x3 = -m*ddx2, -k*x2+2*k*x3-k*x4 = -m*ddx3, -k*x3+2*k*x4-k*x5 = -m*ddx4, -k*x4+2*k*x5-k*x6 = -m*ddx5, -k*x5+2*k*x6-k*x7 = -m*ddx6, -k*x6+2*k*x7-k*x8 = -m*dd...
vector([2*k*x1-k*x2 = -m*ddx1, -k*x1+2*k*x2-k*x3 = -m*ddx2, -k*x2+2*k*x3-k*x4 = -m*ddx3, -k*x3+2*k*x4-k*x5 = -m*ddx4, -k*x4+2*k*x5-k*x6 = -m*ddx5, -k*x5+2*k*x6-k*x7 = -m*ddx6, -k*x6+2*k*x7-k*x8 = -m*dd...
vector([2*k*x1-k*x2 = -m*ddx1, -k*x1+2*k*x2-k*x3 = -m*ddx2, -k*x2+2*k*x3-k*x4 = -m*ddx3, -k*x3+2*k*x4-k*x5 = -m*ddx4, -k*x4+2*k*x5-k*x6 = -m*ddx5, -k*x5+2*k*x6-k*x7 = -m*ddx6, -k*x6+2*k*x7-k*x8 = -m*dd...
vector([2*k*x1-k*x2 = -m*ddx1, -k*x1+2*k*x2-k*x3 = -m*ddx2, -k*x2+2*k*x3-k*x4 = -m*ddx3, -k*x3+2*k*x4-k*x5 = -m*ddx4, -k*x4+2*k*x5-k*x6 = -m*ddx5, -k*x5+2*k*x6-k*x7 = -m*ddx6, -k*x6+2*k*x7-k*x8 = -m*dd...

>    K22:=genmatrix(convert(de22,list),[seq(x||jj,jj=1..22)]);

K22 := matrix([[2*k, -k, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [-k, 2*k, -k, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, -k, 2*k, -k, 0, 0, 0, 0, 0, 0, 0, 0, 0,...

>    K22:=subs(k=1.25e8,evalm(K22));

K22 := matrix([[.250e9, -.125e9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [-.125e9, .250e9, -.125e9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, -.125e9, .250e9, -...

K22は22×22の行列になります。
 
質量行列は, 前のようにmakesysでも作れますが, こ
こでは簡単に帯行列として定義してしまいます。

>    M22:=band([m],22);

M22 := matrix([[m, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, m, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, m, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...

>    M22:=subs(m=5000,evalm(M22));

M22 := matrix([[5000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 5000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 5000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...

こちらもMATLABリンクで固有周波数とモードベクトルを求めます。

>    setvar("K22",K22);

>    setvar("M22",M22);

>    evalM("[P22,W22]=eig(K22,M22)");

>    P22:=getvar("P22"):

>    W22:=getvar("W22"):

>    P22:=convert(P22,array):

>    for jj from 1 to 22 do
  v||jj:=col(P22, jj): w||jj:=W22[jj,jj]:
end do:

>    Matlab[closelink]();  # MATLABリンクを終了

最初の例題と同様, モードベクトルをとりだし, プロットやアニメーションを行うことができます。さまざまなモードベクトルをご覧下さい。

>    listplot(v5,view=[0..22,-0.01..0.01],title=cat("w=",convert(w5,string)),
style=point,symbol=circle);

[Maple Plot]

>    listplot(v3,view=[0..22,-0.01..0.01],title=cat("w=",convert(w3,string)),
style=point,symbol=circle);

[Maple Plot]

最後に V5 のモードベクトルを使って, アニメーションを作ってみましょう。

>    AA:=convert(evalm(v5*sin(2*Pi/10*theta)),list):

>    seq(assign(AAA[theta]=listplot(AA,style=point,symbol=circle,axes=boxed)),theta=1..10);

>    plots[display](seq(AAA[i],i=1..10),insequence=true,view=[0..22,-0.01..0.01]);

[Maple Plot]

  出力されたアニメーション画面をクリックすると, 画面の上部にあるツールバーで再生を行うことができます。   (連続再生)ボタンをクリックしてから  (再生) ボタンをクリックすると, 振動を連続的に見ることができます。

以上のように, Maple と MATLAB で数値を受け渡すことにより, 数式モデルの作成と整理 → 係数行列の生成 → MATLABでの数値計算 という作業がスムーズに行えるようになります。