mat命令と複素数計算

 投稿日：2008年11月21日(金)15時11分26秒
  大熊　です。

前回　NO95 NOの御回答有難うございました。山中さんを中山さんと間違えて投稿しました。失礼を御許しください。
その後だいぶ10進BASICを進めてますが、下記の不具合でストップしてます。
MAT T=INV(A)　が出来ないのです。


OPTION ARITHMETIC COMPLEX
LET j=SQR(-1)
LET R1=5000
LET R2=5000
LET C1=0.1*10^( -6 )
LET C2=0.1*10^( -6 )
LET F=100
LET SZ=0.001
LET NP=4
LET ZC1=1/(2*PI*F*C1)
LET ZC2=1/(2*PI*F*C2)
OPTION BASE 1
PRINT "ZCI=";ZC1

DIM A(NP,NP),B(NP,NP),T(NP,NP),E1(NP)

LET A(1,2)=1/R1+SZ*j
LET A(2,3)=1/R2+SZ*j
LET A(2,4)=SZ+1/ZC1*j
LET A(3,4)=SZ+1/ZC2*j
MAT PRINT A

LET E1(1)=1
LET E1(4)=0
PRINT"下記のごとくMAT PRINT E1とやると横に一文字になる。"
MAT PRINT E1
PRINT"下記のごとくMAT PRINT USING REPEAT$RIで縦に一文字並び良

好。"
MAT PRINT USING REPEAT$(" #.#### ",1):E1
MAT T=INV(A)


STOP

早速ですが、上記の文を作り、マトリクス[A]を作り実行すると、最後の　MAT T=INV(A)で
「EXTYPE 3009 引数が定義外の値」となり、ストップします。
どこが不良の原因でしょうか。
また、一列のE(5)などを作り、PRINT E 　をやると、横一列に表示します。
四端子網の　A*E 等は大丈夫でしょうか。縦に変更してT=TRAN(E) そしてA*T でしょうか。 		 

Re: mat命令と複素数計算

 投稿者：山中和義  投稿日：2008年11月21日(金)16時36分32秒
  > No.104[元記事へ]

大熊　正さんへのお返事です。

> 早速ですが、上記の文を作り、マトリクス[A]を作り実行すると、最後の　MAT T=INV(A)で
> 「EXTYPE 3009 引数が定義外の値」となり、ストップします。
> どこが不良の原因でしょうか。

行列Aが、逆行列を持たない行列だからです。
Aは筆算で逆行列は存在しますか？


> また、一列のE(5)などを作り、PRINT E 　をやると、横一列に表示します。
> 四端子網の　A*E 等は大丈夫でしょうか。縦に変更してT=TRAN(E) そしてA*T でしょうか。

DIM A(3,3)
DIM X(3),B(3)
MAT B=A*X !(3行,3列)(3行,1列)=(3行,1列)として計算される
MAT PRINT B !横へ
MAT B=X*A !(1行,3列)(3行,3列)=(1行,3列)として計算される
MAT PRINT B !横へ
END
この場合、X,Bはベクトル扱いになります。
X,Bを行列として扱う場合は、（X,Bに対してTRNを適用する場合）
行または列のみの行列は、たとえば
　３行１列なら DIM B(3,1)
　１行３列なら DIM B(1,3)
としてください。 		 

Re: mat命令と複素数計算

 投稿者：SECOND  投稿日：2008年11月23日(日)07時44分19秒
  > No.104[元記事へ]

大熊　正さんへのお返事です。

> MAT T=INV(A)　が出来ないのです。

IF DET(A)<>0 THEN !　行列式｜A｜の値
   MAT T=INV(A)
ELSE
   PRINT "A は、逆行列を持たない"
END IF


※十進BASIC の、１次元配列、行ベクトルと列ベクトル

行列 (a)
┌              ┐
│ 1.000   2.000│
│ 3.000   4.000│
└              ┘
ベクトル (v1)
(  1.000   2.000 )　･･･この状態は、行か、列かが、不定になっている。

MAT v2=a*v1　･･･右へ書けば、列ベクトル v1 として計算される。
┌              ┐┌      ┐
│ 1.000   2.000││ 1.000│
│ 3.000   4.000││ 2.000│
└              ┘└      ┘
v2＝(  5.000  11.000 )

MAT v2=v1*a　･･･左へ書けば、行ベクトル v1 として計算される。
┌              ┐┌              ┐
│ 1.000   2.000││ 1.000   2.000│
└              ┘│ 3.000   4.000│
                  └              ┘
v2＝(  7.000  10.000 ) 		 

Re: mat命令と複素数計算

 投稿者：島村1243  投稿日：2008年11月23日(日)10時38分16秒
  > No.104[元記事へ]

大熊　正さんへのお返事です。

> 大熊　です。
> ***中略***
> その後だいぶ10進BASICを進めてますが、下記の不具合でストップしてます。
> MAT T=INV(A)　が出来ないのです。
> ***中略***
> LET A(1,2)=1/R1+SZ*j
> LET A(2,3)=1/R2+SZ*j
> LET A(2,4)=SZ+1/ZC1*j
> LET A(3,4)=SZ+1/ZC2*j
> MAT PRINT A
> ***中略***
> 早速ですが、上記の文を作り、マトリクス[A]を作り実行すると、最後の　MAT T=INV(A)で
> 「EXTYPE 3009 引数が定義外の値」となり、ストップします。
> どこが不良の原因でしょうか。

アドミタンス行列[A]と節点法を使って、電気回路網の電流[I]=[A][E]、インピーダンス[Z]=INV(A)を計算するプログラムを作成するつもりの様ですね。
プログラム原稿を見ると、下記節点間アドミタンスが未設定であることがエラーの原因です。

未設定の節点間アドミタンス
A(2,1)!=A(1,2)
A(3,2)!=A(2,3)
A(4,3)!=A(3,4)

したがって「INV(A)」よりも前に
A(2,1)=A(1,2)
A(3,2)=A(2,3)
A(4,3)=A(3,4)

を追記すればエラーは出ません。
なお、本来は対地間アドミタンス｛A(1,1),A(2,2),A(3,3),A(4,4)｝や、アドミタンスが接続されない節点間のアドミタンスは０と 設定すべきです（本件ではBASICアプリケーションが未指定の節点間アドミタンス＝０と自動初期値設定をしてしまうのでトラブルにはなりませんが注意が 必要です）。 		 

mat命令と複素数計算

 投稿日：2008年11月24日(月)18時50分22秒
  大熊　です。
山中様、島村様、そしてSECOND様
色いろの御教授　本当に有難うございました。

その総てが、現在の私に理解できたとは、到底思ってませんが、
10進BASICが更に好きになったことは確かです。「持つべきは、
先達なり」・・・と感謝いたしております。

そこで、電気回路{A}では、コンデンサーなどのアドミッタンス
が周波数（Ｆ）特性を持ちます、更に電源[B}も同様です。

（１）マトリクス[A]の中にコンデンサーなどのアドミッタンス
　　　をスマートに入れる方法。
（２）同じく　電源のマトリクス[B}に、SIN（F,T)表示等で
　　　スマートに入れる方法。
（３）電源の周波数と大きさ、出力の関係、将来のグラフ化に
　　　備え、INV(A)*B　等とした上で更に、全体を
　　　　　FOR F=100 TO 10000 STEP 100・・「INB (A)*B」・
　　　　　・・・NEXT F
       等とやるのでしょうか、全体の方法が私にはまだ見え
　　　　てません。
（４）最終的には、片対数表のグラフ表示に成るのでしょうが、
　　　これを　　　やった、参考資料や、文献などがありまし
　　　たら御教えください。
　　　ＣＲ　一段の回路でも現在の私には、大変参考になります。
（５）実は、最終的には「有限要素法」にまでたどり着きた
　　　いのですが、　　　10進BASICでこれをやった、先達の文献
　　　などありましたら御教えください。従来の　[N88 BASIC ??]
　　　でやった参考書はあるのですが、真似してプログラムすると
　　　「文法の相違か、あちこちでつっかえ全く動かず」
　　　現在は、諦めている状況です。

「もっと、自分で苦労し、真面目にやれ」とのお叱りの言葉も
　きこえますが、なるべく早く「初歩の段階」を済ませ先に
　行きたいので　よろしく御願いします。

敬具 		 

Re: mat命令と複素数計算

 投稿者：山中和義  投稿日：2008年11月25日(火)11時29分9秒
  > No.114[元記事へ]

大熊　正さんへのお返事です。

> （１）マトリクス[A]の中にコンデンサーなどのアドミッタンスをスマートに入れる方法。
> （２）同じく　電源のマトリクス[B}に、SIN（F,T)表示等でスマートに入れる方法。

１要素ずつ代入文で設定するのが基本になります。
サブルーチンや関数を使って、最小限の要素を指示することも可能です。（サンプル参照）

RC回路、これでいいのでしょうか？　具体的に提示した方が問題解決が早いと思います。


> （４）最終的には、片対数表のグラフ表示に成るのでしょうが、

周波数特性のグラフですか？


> （５）実は、最終的には「有限要素法」にまでたどり着きたいのですが、

tについての微分方程式を解くということですか？


２端子対回路（４端子回路）として処理したサンプル
1000 !２端子対回路（４端子回路）
1010 !　i1→┌──┐i2→
1020 !　a ─┤A　B├─ b
1030 !E1↑　│　　│　↑E2
1040 !　a'─┤C　D├─ b'
1050 !　　　└──┘
1060 !基本行列Fを用いて
1070 ! (E1)=(A B)(E2)
1080 ! (i1) (C D)(i2)
1090
1100 OPTION ARITHMETIC COMPLEX !複素数を扱う
1110
1120 LET j=SQR(-1) !虚数単位　※電気系はｊを使う
1130
1140 !●交流回路
1150
1160 LET f=60 !周波数[Hz]
1170 DEF w=2*PI*f !角周波数ω
1180 DEF w2f=w/(2*PI) !ωからfを求める
1190
1200 DEF H2Ohm(L)=j*w*L ![H]を[Ω]へ
1210 DEF F2Ohm(C)=1/(j*w*C) ![F]を[Ω]へ
1220 DEF xL(L)=w*L !誘導リアクタンス
1230 DEF xC(C)=1/(w*C) !容量リアクタンス
1240
1250 SUB DispS(z) !複素数をS表示する　※スタインメッツ（Steinmetz）
1260    PRINT ABS(z);
1270    IF ABS(z)<>0 THEN
1280       IF arg(z)<>0 THEN PRINT "∠";DEG(arg(z));"°";
1290    END IF
1300    PRINT
1310 END SUB
1320
1330 FUNCTION S2COMPLEX(l,th) !S表示（極座標形式）を複素数へ
1340    LET S2COMPLEX=COMPLEX(l*COS(RAD(th)),l*SIN(RAD(th)))
1350 END FUNCTION
1360 FUNCTION i2COMPLEX(im,th) !瞬時値式を複素数へ　※最大値、初期位相
1370    LET i2COMPLEX=S2COMPLEX(im/SQR(2),th) !実効値、初期位相
1380 END FUNCTION
1390
1400
1410 !●基本２端子対回路
1420 !　─Z─　の場合　F=(1 Z)
1430 !　─-─　　　　　　(0 1)
1440 SUB seriesF(Z,F(,)) !直列挿入
1450    MAT F=IDN
1460    LET F(1,2)=Z !インピーダンス
1470 END SUB
1480 !　─┬─　の場合　F=(1　 0)
1490 !　　Z 　　　　　　　(1/Z 1)
1500 !　─┴─
1510 SUB shuntF(Z,F(,)) !並列挿入
1520    WHEN EXCEPTION IN
1530       MAT F=IDN
1540       LET F(2,1)=1/Z !アドミタンス
1550    USE
1560       PRINT "０では割れません。"
1570       STOP
1580    END WHEN
1590 END SUB
1600
1610 !●パラメータの相互変換　※一部
1620 SUB F2Z(F(,), Z(,)) !ＦパラメータをＺパラメータへ
1630    LET Z(1,1)=F(1,1) !A
1640    LET Z(1,2)=DET(F) !AD-BC
1650    LET Z(2,1)=1
1660    LET Z(2,2)=F(2,2) !D
1670    WHEN EXCEPTION IN
1680       MAT Z=(1/F(2,1))*Z !(1/C)倍
1690    USE
1700       PRINT "Ｚパラメータは存在しません。"
1710       STOP
1720    END WHEN
1730 END SUB
1740 SUB F2Y(F(,),Y(,)) !ＦパラメータをＹパラメータへ
1750    LET Y(1,1)=F(2,2) !D
1760    LET Y(1,2)=-DET(F) !BC-AD
1770    LET Y(2,1)=-1
1780    LET Y(2,2)=F(1,1) !A
1790    WHEN EXCEPTION IN
1800       MAT Y=(1/F(1,2))*Y !(1/B)倍
1810    USE
1820       PRINT "Ｙパラメータは存在しません。"
1830       STOP
1840    END WHEN
1850 END SUB
1860 !-------------------- ここまでがサブルーチン
1870
1880
1890 DIM mF(2,2),mZ(2,2),mY(2,2) !F,Z,Yパラメータ
1900 DIM vi(2),vo(2) !電流や電圧のベクトル
1910 DIM T1(2,2),T2(2,2) !作業用
1920
1930
1940 !●１次ＲＣ回路、直列ＲＣ回路、ローパスフィルタ
1950 !i1→　　　i2→
1960 !　─R─┬─
1970 !E1↑　 C　↑E2
1980 !　─-─┴─
1990 !i1=3*SQR(2)*SIN(377*t)[A]、R=30[Ω]、C=66.3[μF]
2000
2010 LET R=30
2020 LET C=66.3e-6
2030
2040
2050 !●Ｆパラメータ、基本行列
2060 MAT mF=IDN !単位行列
2070
2080 CALL seriesF(R,T1)
2090 MAT mF=mF*T1 !縦続接続
2100
2110 CALL shuntF(F2Ohm(C),T2)
2120 MAT mF=mF*T2
2130
2140 PRINT "Ｆパラメータ"
2150 MAT PRINT mF;
2160
2170
2180 !●Ｚパラメータ、インピーダンス行列 V=Z*I
2190 CALL F2Z(mF,mZ)
2200 PRINT "Ｚパラメータ"
2210 MAT PRINT mZ;
2220
2230
2240 !●(E1)=[F](E2)より
2250 !　(i1)　　(i2)
2260 LET vi(2)=i2COMPLEX(3*SQR(2),0) !3[A}
2270 PRINT "i1=";
2280 CALL DispS(vi(2))
2290 LET vi(1)=mZ(1,1)*vi(2) !E=R*i、150∠-53.1°[V]
2300 PRINT "E1=";
2310 CALL DispS(vi(1))
2320 PRINT
2330
2340 MAT T1=INV(mF) !出力側を算出する
2350 MAT vo=T1*vi
2360
2370 PRINT "E2=";
2380 CALL DispS(vo(1))
2390 PRINT "i2=";
2400 CALL DispS(vo(2))
2410 PRINT
2420
2430
2440 !●Ｙパラメータ、アドミタンス行列 I=Y*V
2450 CALL F2Y(mF,mY)
2460 PRINT "Ｙパラメータ"
2470 MAT PRINT mY;
2480
2490
2500 END
 

Re: mat命令と複素数計算

 投稿者：山中和義  投稿日：2008年11月26日(水)08時19分38秒
  > No.114[元記事へ]

大熊　正さんへのお返事です。
> （１）マトリクス[A]の中にコンデンサーなどのアドミッタンスをスマートに入れる方法。
> （２）同じく　電源のマトリクス[B}に、SIN（F,T)表示等でスマートに入れる方法。

節点電位法によるサンプルです。
以前第１掲示板に投稿した直流回路の改修版になります。
1000 !電気回路シミュレーション（節点電位法）　交流回路
1010
1020 !・各節点の電位を表示する
1030 !・各素子への電流、電位を表示する
1040
1050 OPTION ARITHMETIC COMPLEX
1060
1070 LET j=SQR(-1) !虚数単位　※電気系はｊを使う
1080
1090 LET f=60 !周波数[Hz]
1100 DEF w=2*PI*f !角周波数ω
1110
1120 DEF H2Ohm(L)=j*w*L ![H]を[Ω]へ
1130 DEF F2Ohm(C)=1/(j*w*C) ![F]を[Ω]へ
1140 DEF xL(L)=w*L !誘導リアクタンス
1150 DEF xC(C)=1/(w*C) !容量リアクタンス
1160
1170 SUB DispS(z) !複素数をS表示する　※スタインメッツ（Steinmetz）
1180    PRINT ABS(z);
1190    IF ABS(z)<>0 THEN
1200       IF arg(z)<>0 THEN PRINT "∠";DEG(arg(z));"°";
1210    END IF
1220    !PRINT
1230 END SUB
1240
1250 FUNCTION S2COMPLEX(l,th) !S表示（極座標形式）を複素数へ
1260    LET S2COMPLEX=COMPLEX(l*COS(RAD(th)),l*SIN(RAD(th)))
1270 END FUNCTION
1280 FUNCTION i2COMPLEX(im,th) !瞬時値式を複素数へ　※最大値、初期位相
1290    LET i2COMPLEX=S2COMPLEX(im/SQR(2),th) !実効値、初期位相
1300 END FUNCTION
1310 !-------------------- ここまでがサブルーチン
1320
1330 !---------- ↓↓↓↓↓ ----------
1340 LET Ns=0 !電圧源の数
1350 LET Nd=3 !節点の数
1360 !---------- ↑↑↑↑↑ ----------
1370
1380 LET N=Nd+Ns
1390
1400 !●キルヒホッフの電流則より、節点方程式を組み立てる
1410 DIM A(N,N),x(N),b(N) !連立方程式 Ax=b
1420 MAT A=ZER
1430 MAT b=ZER
1440
1450
1460 DIM el_tmp$(100),ev_tmp(100),nd1_tmp(100),nd2_tmp(100) !素子の属性
1470 LET K=0
1480
1490 SUB AddElements(el$,ev,nd1,nd2) !回路を記録する
1500    LET K=K+1 !連番で回路を記録する　※注意. 方程式の成分番号と一致しない
1510    LET el_tmp$(K)=el$
1520    LET ev_tmp(K)=ev
1530    LET nd1_tmp(K)=nd1
1540    LET nd2_tmp(K)=nd2
1550
1560
1570    !連立方程式を組み立てる
1580    ! ┌　 │　 ┐┌　┐ ┌　┐
1590    ! │Ｇ │±1││V │=│I │
1600    !───┼─────────
1610    ! │±1│-Zp││Ip│ │Ep│
1620    ! └　 │　 ┘└　┘ └　┘
1630
1640    SELECT CASE UCASE$(el$(1:1)) !素子に応じて
1650    CASE "V" !電圧源なら
1660       LET p=VAL(el$(2:LEN(el$)))+Nd !番号を得る
1670       LET A(nd1,p)=A(nd1,p)-1 !電流Ipが節点iから節点jへ流れたとして、(Vi-Vj)-Ip*Zp=Ep
1680       LET A(p,nd1)=A(p,nd1)-1
1690       LET A(nd2,p)=A(nd2,p)+1
1700       LET A(p,nd2)=A(p,nd2)+1
1710       LET A(p,p)=A(p,p)+0 !内部抵抗Zpは０とする
1720       LET b(p)=b(p)+ev !起電力
1730
1740    CASE "I" !電流源なら
1750       LET b(nd1)=b(nd1)-ev
1760       LET b(nd2)=b(nd2)+ev
1770
1780    CASE ELSE !素子なら
1790       LET Gij=1/ev
1800       !対角成分　※節点に接続された素子（アドミタンス）の和
1810       LET A(nd1,nd1)=A(nd1,nd1)+Gij
1820       LET A(nd2,nd2)=A(nd2,nd2)+Gij
1830       !その他の成分　※節点に接続された素子（アドミタンス）に-1をかけたものの和
1840       LET A(nd1,nd2)=A(nd1,nd2)-Gij
1850       LET A(nd2,nd1)=A(nd2,nd1)-Gij
1860
1870    END SELECT
1880 END SUB
1890 !-------------------- ここまでがサブルーチン
1900
1910 !---------- ↓↓↓↓↓ ----------
1920
1930 !●回路図　※１次ＲＣ回路、直列ＲＣ回路、ローパスフィルタ
1940 !　─２─R1─３─
1950 !　　│　　　│
1960 !　　i3　　　C2
1970 !　　│　　　│
1980 !　─１───┴─
1990 !　　│
2000 !　　≡アース
2010 !　i3=3*SQR(2)*SIN(377*t)[A]、R1=30[Ω]、C2=66.3[μF]
2020
2030 !素子: Rn,Vn,In、n:番号（連番）　※２文字目以降は番号
2040 !値:
2050 !端子番号（起点）: １以上の値　※節点
2060 !端子番号: １以上の値　※節点
2070
2080 CALL AddElements("R1",30,2,3) !30[Ω]、枝路電流は2→3と仮定する
2090 CALL AddElements("C2",F2Ohm(66.3e-6),3,1) !C=66.3[μF]
2100 CALL AddElements("i3",i2COMPLEX(3*SQR(2),0),1,2) !3[A]
2110
2120 !※電圧源の番号は１からの連番　例. CALL AddElements("V1",?,?,?) !?[V]
2130 !なし
2140
2150 LET GND=1 !※アース
2160
2170 !---------- ↑↑↑↑↑ ----------
2180
2190
2200 FOR i=1 TO Nd !結線されていない節点　1*Vi=0
2210    IF A(i,i)=0 THEN LET A(i,i)=1
2220 NEXT i
2230 LET A(GND,GND)=0 !電位を０とする
2240
2250 MAT PRINT A;
2260 MAT PRINT b;
2270
2280
2290 DIM Ai(N,N) !連立方程式を解く
2300 MAT Ai=INV(A)
2310 MAT x=Ai*b
2320
2330
2340 FOR i=1 TO Nd !各節点の電位、流れ込む電流を表示する
2350    PRINT "節点";STR$(i);":";
2360    CALL DispS(x(i))
2370    PRINT "[V] ,";
2380    CALL DispS(b(i))
2390    PRINT "[A]"
2400 NEXT i
2410 PRINT
2420
2430 FOR i=1 TO K-Ns !各素子に流れる電流、電位を表示する
2440    PRINT el_tmp$(i);":";
2450    LET t$=el_tmp$(i)(1:1)
2460    IF UCASE$(t$)="I" THEN !電流源なら
2470       CALL DispS(ev_tmp(i))
2480       PRINT "[A]",
2490       CALL DispS(x(nd2_tmp(i))-x(nd1_tmp(i)))
2500       PRINT "[V]"
2510    ELSE
2520       CALL DispS((x(nd1_tmp(i))-x(nd2_tmp(i)))/ev_tmp(i))
2530       PRINT "[A]",
2540       CALL DispS(x(nd1_tmp(i))-x(nd2_tmp(i)))
2550       PRINT "[V]"
2560    END IF
2570 NEXT i
2580 PRINT
2590
2600 FOR i=1 TO Ns !電圧源に流れる電流、電位を表示する
2610    PRINT el_tmp$(K-Ns+i);":";
2620    CALL DispS(x(Nd+i))
2630    PRINT "[A]",
2640    CALL DispS(ev_tmp(K-Ns+i))
2650    PRINT "[V]"
2660 NEXT i
2670
2680
2690 END
 

戻る