Chapter 版 Maxima を用いた LC のインピーダンス測定について [ 目的 ] 電気通信大学 先進理工学科の2 年次後期に実施される電気 電子回路実験において L,C のインピーダンス測定を実施している この実験項目について 無料ソフトの Maxima を用い

Chapter 2 [目的] 電気通信大学・先進理工学科の２年次後期に実施される電気・電子回路実験において L,C のインピーダンス測定を実施している。この実験項目について、無料ソフトの Maxima を用い ることで、理論解析と実験値の比較が可能である。また、近年のパソコンの性能の向上により、 Maxima の実行処理速度が大幅に改善された。 Maxima を用いて、計算方法と計算結果を示すことで、実験レポートの考察のヒントにして、実 験内容の理解を深めることが目的である。 [目次] 以下の8 節から構成されています。

2-1

コイルにおけるインピーダンスの大きさと周波数の関係

2-2

実験データとコイルにおけるインピーダンスと周波数の関係

2-3

コイルにおけるインピーダンスの位相と周波数の関係

2-4

実験データとコイルにおけるインピーダンスの位相と周波数の関係

2-5

コンデンサにおけるインピーダンスの大きさと周波数の関係

2-6

実験データとコンデンサにおけるインピーダンスの大きさと周波数の関係

2-7

コンデンサにおけるインピーダンスの位相と周波数の関係

2-8

実験データとコンデンサにおけるインピーダンスの位相と周波数の関係 課題について 課題解答例 1-A)コイルのインピーダンスの大きさにおける誤差の評価 1-B)コイルの位相計算に含まれる誤差の評価について

Maxima を用いた LC のインピーダンス測定について

2016.10.14 版

2-1 コイルにおけるインピーダンスの大きさと周波数の関係

[目的] 周波数依存性のグラフを表示する。

[結果]

変数をすべて消去する。

(%i1)

kill(all);

直列抵抗:r が存在するコイルの

インピーダンスを z とする。

(%i1)

z:r+%i*w*L;

変数は、正の値と仮定する。

(%i2)

assume(r>0,L>0);

Cabs 関数をもちいて、大きさを計算する。

(%i3)

absL:cabs(z);

内部抵抗：r=0.1ohm、インダクタンス：L=1mH を代入する。

(%i4)

ev(absL,r:0.1,L:0.001);

角周波数：w を周波数：f に変換する。

(%i5)

ratsubst(2*%pi*f,w,%o4);

Fig.1 計算する回路図

周波数の関数として、インピーダンスを定義する。

(%i6)

define(AzL(f),%o5);

以上で、周波数依存する関数が定義できたので、

インピーダンスと周波数の関係のグラフを表示する。

グラフの軸を log-log で表示する。

(%i7)

wxplot2d(AzL(f),[f,100,50*10^3],

[logx],[logy],

[xlabel,"Frq"],

[ylabel,"AzL"],

[gnuplot_preamble,"set grid"]);

Log-log グラフで表示すると、

直線状のグラフが確認できます。

Fig.2 計算結果のグラフ

2-2 実験データとコイルにおけるインピーダンスと周波数の関係 [目的] 理論値、実験結果と LCR メータの値を比較する。 [結果] 変数を初期化する。 (%i1) kill(all); インピーダンスを定義する。 (%i1) z:r+%i*w*L; 変数を正値と定義する。 (%i2) assume(r>0,L>0); 複素数の大きさを求める関数を使用 (%i3) absL:cabs(z); 値を代入 (%i4) ev(absL,r:0.1,L:0.001); 周波数に変換 (%i5) ratsubst(2*%pi*f,w,%o4);

Fig.1 計算する回路図

周波数の関数を定義 (%i6) define(AzL(f),%o5); 理論値を表示する。 (%i7) wxplot2d(AzL(f),[f,100,50*10^3], [logx],[logy], [xlabel,"Frq"], [ylabel,"AzL"], [gnuplot_preamble,"set grid"]); 実験データをリスト型で入力する。 (%i8) datax:[500, 1000, 2000, 5000, 10000, 20000, 50000, 100000]; (%i9) datay :[ 3.28, 6.20, 12.04, 30.22, 60.12, 121.63, 305.34, 709.08];

Fig.2 理論グラフ

実験データを表示する。 (%i10) wxplot2d([discrete,datax,datay],[x,100,100*10^3], [logx],[logy],[style,points],[color,red],[gnuplot_preamble,"set grid"]); 理論値と実験値のグラフを表示する。 (%i11) wxplot2d([AzL(x),[discrete,datax,datay]], [x,100,100*10^3], [logx],[logy], [xlabel, "Frequency[Hz]"], [ylabel,"Phase[Deg.]"], [legend,false], [color,blue,red], [style,lines,points], [point_type,circule], [gnuplot_preamble,"set grid"]);

Fig.3 実験値グラフ

Fig.4 理論と実験値グラフ

LCR メーター値をリスト型として入力する。 (%i12) datax1:[500, 1000, 2000, 5000, 10000, 20000, 50000, 100000]; (%i13) datay1 :[ 3.1765, 6.1353, 12.127, 30.069, 59.868, 119.332, 297.668, 596.584]; 確認のためグラフで表示する。 (%i14) wxplot2d([discrete,datax1,datay1],[x,100,100*10^3], [logx],[logy],[style,points],[color,green],[gnuplot_preamble,"set grid"]); 最後に、実験データと理論値と LCR メーター値を同時に表示する。 (%i15) wxplot2d([AzL(x),[discrete,datax,datay],[discrete,datax1,datay1]], [x,100,100*10^3], [logx],[logy], [xlabel, "Frequency[Hz]"], [ylabel,"Phase[Deg.]"], [legend,false], [color,blue,red,green], [style,lines,points,points], [point_type,circule],

Fig.5 LCR 値のグラフ

Fig.6 理論、LCR 値、実験値グラフ

[Maxima ファイルの加工について]

以上の結果は、Maxima ver.5.23 を使用しました。

この version での日本語入力について、

2-3 コイルにおけるインピーダンスの位相と周波数の関係

[目的] 位相の周波数依存性グラフを表示する。

[結果]

すべての変数を初期化する。

(%i1)

kill(all);

インピーダンスｚを定義する。

(%i1)

z:r+%i*w*L;

変数を正の値と定義する。

(%i2)

assume(r>0,L>0);

位相を計算する関数 carg 関数を用いる。

(%i3)

thL:carg(z);

内部抵抗：r=0.1ohm、インダクタンス：L=1mH を代入する。

(%i4)

ev(thL,r:0.1,L:0.001);

角周波数：w を周波数：ｆに変換する。

(%i5)

ratsubst(2*%pi*f,w,%o4);

Fig.1 計算する回路図

周波数ｆの関数として、位相を定義する。

(%i6)

define(thetaL(f),%o5);

以上までに、関数の定義ができたので、

位相と周波数の関係のグラフを表示する。

(%i7)

wxplot2d(thetaL(f)*180/%pi,[f,100,50*10^3],

[logx],

[xlabel,"Frq"],

[ylabel,"AzL"],

[gnuplot_preamble,"set grid"]);

低周波で、純抵抗の性質に近似できるので、

位相がゼロに近くなる。

Fig.2 計算結果のグラフ

2-4 実験データとコイルにおけるインピーダンスの位相と周波数の関係

[目的] 理論値、実験結果と LCR メータの結果を比較する。 [結果] 変数の初期化 (%i1) kill(all); インピーダンスを定義 (%i1) z:r+%i*w*L; 変数は、正の値とする。 (%i2) assume(r>0,L>0); 複素数の位相を求める関数を使用する。 (%i3) thL:carg(z); 変数に、値を代入する。 内部抵抗：r=0.1ohm、インダクタンス：L=1mH を代入する。 (%i4) ev(thL,r:1.0,L:0.001); 変数を周波数に変換 (%i5) ratsubst(2*%pi*f,w,%o4); 周波数の関数を定義する。 (%i6) define(thetaL(f),%o5);

Fig.1 計算する回路図

理論線を表示する。 (%i7) wxplot2d(thetaL(f)*180/%pi,[f,100,50*10^3], [logx], [xlabel,"Frq"], [ylabel,"AzL"], [gnuplot_preamble,"set grid"]); 実験データをリストで定義して入力する。 Data from 1313068,sagara

(%i8) datax:[500, 1000, 2000, 5000, 10000, 20000, 50000, 100000];

(%i9) datay :[ 75.6, 82.8, 86.4,90.0, 86.4, 89.3, 93.6, 90.0];

実験値を表示する。 (%i10) wxplot2d([discrete,datax,datay],[x,100,100*10^3], [logx],[style,points],[color,red],[gnuplot_preamble,"set grid"]); 理論値と実験値を表示する。 (%i11) wxplot2d([thetaL(x)*180/%pi,[discrete,datax,datay]], [x,100,100*10^3], [logx], [xlabel, "Frequency[Hz]"], [ylabel,"Phase[Deg.]"], [legend,false], [color,blue,red], [style,lines,points], [point_type,circule], [gnuplot_preamble,"set grid"]);

Fig.3 実験値のグラフ

Fig.4 理論値と実験値のグラフ

LCR メーターの測定結果をリスト型で入力する。 LCR-meter Data from 1313068,sagara

(%i12) datax1:[500, 1000, 2000, 5000, 10000, 20000, 50000, 100000]; (%i13) datay1 :[ 72.838, 81.024, 85.237,87.767, 88.748, 89.241, 89.462, 89.395]; 実験値、理論値と LCR メーター値を同時に表示する。 (%i14) wxplot2d([thetaL(x)*180/%pi,[discrete,datax,datay],[discrete,datax1,datay1]], [x,100,100*10^3], [logx], [xlabel, "Frequency[Hz]"], [ylabel,"Phase[Deg.]"], [legend,false], [color,blue,red,green], [style,lines,points,points], [point_type,circule], [gnuplot_preamble,"set grid"]);

Fig.5 すべての結果のグラフ

2-5 コンデンサにおけるインピーダンスの大きさと周波数の関係

[目的]

周波数依存性のグラフを表示する。

[結果]

変数をすべて消去する。

(%i1)

kill(all);

直列抵抗:r が存在するコンデンサの

インピーダンスを z とする。

(%i1)

z:r+1/(%i*w*C);

変数は、正の値と仮定する。

(%i2)

assume(r>0,C>0);

Cabs 関数をもちいて、大きさを計算する。

(%i3)

absC:cabs(z);

内部抵抗：r=0.01ohm、キャパシタンス：C=1uF の値を代入する。

(%i4)

ev(absC,r:0.01,C:0.1*10^-6);

角周波数：w を周波数：f に変換する。

(%i5)

ratsubst(2*%pi*f,w,%o4);

Fig.1 計算する回路図

周波数の関数として、インピーダンスを定義する。

(%i6)

define(AbsC(f),%o5);

以上で、周波数依存する関数が定義できたので、

インピーダンスと周波数の関係のグラフを表示する。

(%i7)

wxplot2d(AbsC(f),[f,100,50*10^3],

[logx],[logy],

[xlabel,"Frq"],

[ylabel,"AzL"],

[gnuplot_preamble,"set grid"]);

周波数の増加により、インピーダンスが減少することが確認できる。

Fig.2 計算結果のグラフ

2-6 実験データとコンデンサにおけるインピーダンスの大きさと周波数の関係

[目的] 理論値、実験結果と LCR メータの値を比較する。 [結果] 変数を初期化する。 (%i15) kill(all); インピーダンスを定義する。 (%i1) z:r+1/(%i*w*C); 変数を正値と定義する。 (%i2) assume(r>0,C>0); 複素数の大きさを求める関数を使用 (%i3) absC:cabs(z); 容量 C=0.1uF と等価抵抗 r=0.01ohm の値を代入 (%i4) ev(absC,r:0.01,C:0.1*10^-6); 周波数に変換 (%i5) ratsubst(2*%pi*f,w,%o4);

Fig.1 計算する回路図

周波数の関数を定義 (%i6) define(AbsC(f),%o5); 理論値を表示する。 (%i7) wxplot2d(AbsC(f),[f,100,50*10^3], [logx],[logy], [xlabel,"Frq"], [ylabel,"AzL"], [gnuplot_preamble,"set grid"]); 実験データをリスト型で入力する。 (%i8) datax:[500, 1000, 2000, 5000, 10000, 20000, 50000, 100000]; (%i9) datay :[ 3173.88, 1586.94, 780.02, 316.65, 156.94, 79.44, 31.48, 15.85];

Fig.2 理論グラフ

実験データを表示する。 (%i10) wxplot2d([discrete,datax,datay],[x,100,100*10^3], [logx],[logy],[style,points],[color,red],[gnuplot_preamble,"set grid"]); 理論値と実験値のグラフを表示する。 (%i11) wxplot2d([AbsC(x),[discrete,datax,datay]], [x,100,100*10^3], [logx],[logy], [xlabel, "Frequency[Hz]"], [ylabel,"|Z|[ohm]"], [legend,false], [color,blue,red], [style,lines,points], [point_type,circule], [gnuplot_preamble,"set grid"]);

Fig.3 実験値グラフ

LCR メーター値をリスト型として入力する。 (%i12) datax1:[500, 1000, 2000, 5000, 10000, 20000, 50000, 100000]; (%i13) datay1 :[ 3124.914, 1564.86, 784.146, 314.779, 157.940, 79.303, 31.930, 16.019]; 最後に、実験データと理論値と LCR メーター値を同時に表示する。 (%i14) wxplot2d([AbsC(x),[discrete,datax,datay],[discrete,datax1,datay1]], [x,100,100*10^3], [logx],[logy], [xlabel, "Frequency[Hz]"], [ylabel,"|Z|[ohm]"], [legend,false], [color,blue,red,green], [style,lines,points,points], [point_type,circule], [gnuplot_preamble,"set grid"]);

Fig.5 理論、LCR 値、実験値グラフ

2-7 コンデンサにおけるインピーダンスの位相と周波数の関係

[目的]

周波数依存性のグラフを表示する。

[結果]

すべての変数を初期化する。

(%i1)

kill(all);

インピーダンスｚを定義する。

(%i1)

z:r+1/(%i*w*C);

変数を正の値と定義する。

(%i2)

assume(r>0,C>0);

位相を計算する関数 carg 関数を用いる。

(%i3)

thC:carg(z);

内部抵抗：r=0.01ohm、キャパシタンス：C=0.1uF の数値を代入する。

(%i4)

ev(thC,r:0.01,C:0.1*10^-6);

角周波数：w を周波数：ｆに変換する。

(%i5)

ratsubst(2*%pi*f,w,%o4);

Fig.1 計算する回路図

周波数ｆの関数として、位相を定義する。

(%i6)

define(ThetaC(f),%o5);

以上までに、関数の定義ができたので、

位相と周波数の関係のグラフを表示する。

(%i7)

wxplot2d(ThetaC(f)*180/%pi,[f,100,50*10^3],

[logx],

[xlabel,"Frq"],

[ylabel,"AzL"],

[gnuplot_preamble,"set grid"]);

周波数の増加により、純抵抗に近似できるので、ゼロに近づくことがわかる。

Fig.2 計算結果のグラフ

2-8 実験データとコンデンサにおけるインピーダンスの位相と周波数の関係

[目的] 理論値、実験結果と LCR メータの値を比較する。 [結果] 変数を初期化する。 (%i15) kill(all); インピーダンスを定義する。 (%i1) z:r+1/(%i*w*C); 変数を正値と定義する。 (%i2) assume(r>0,C>0); 複素数の大きさを求める関数を使用 (%i3) thC:carg(z); 容量 C=0.1uF,直列抵抗 r=0.3 ohm の値を代入する。 (%i4) ev(thC,r:0.3,C:0.1*10^-6); 周波数に変換 (%i5) ratsubst(2*%pi*f,w,%o4);

Fig.1 計算する回路図

周波数の関数を定義 (%i6) define(ThetaC(f),%o5); 理論値を表示する。 (%i7) wxplot2d(ThetaC(f)*180/%pi,[f,100,50*10^3], [logx], [xlabel,"Frq"], [ylabel,"Phase[Deg.]"], [gnuplot_preamble,"set grid"]); 実験データをリスト型で入力する。 (%i8) datax:[500,1000,2000,5000,10*10^3,20*10^3,50*10^3,100*10^3]; (%i9) datay:[-86.4,-82.8,-89.3,-86.4,-90,-89.3,-86.4,-86.4];

Fig.2 理論グラフ

実験データを表示する。 (%i10) wxplot2d([discrete,datax,datay],[x,100,1*10^6], [logx],[style,points],[color,red],[gnuplot_preamble,"set grid"]); 理論値と実験値のグラフを表示する。 (%i11) wxplot2d([ThetaC(x)*180/%pi,[discrete,datax,datay]], [x,100,1*10^6], [logx], [xlabel, "Frequency[Hz]"], [ylabel,"Phase[Deg.]"], [legend,false], [color,blue,red], [style,lines,points], [point_type,circule], [gnuplot_preamble,"set grid"]);

Fig.4 理論と実験値グラフ

Fig.3 実験値グラフ

LCR メーター値をリスト型として入力する。 (%i12) datax1:[500,1000,2000,5000,10*10^3,20*10^3,50*10^3,100*10^3]; (%i13) datay1:[-89.813,-89.752,-89.673,-89.545,-89.420,-89.242,-88.881,-88.394]; 最後に、実験データと理論値と LCR メーター値を同時に表示する。 (%i14) wxplot2d([ThetaC(x)*180/%pi,[discrete,datax,datay],[discrete,datax1,datay1]], [x,100,1*10^6], [logx], [xlabel, "Frequency[Hz]"], [ylabel,"Phase[Deg.]"], [legend,false], [color,blue,red,green], [style,lines,points,points], [point_type,circule], [gnuplot_preamble,"set grid"]);

Fig.5 理論、LCR 値、実験値グラフ

Chapter2 課題について 理論との相違を考慮するために、誤差の計算をする。 誤差の計算により、理論値と一致するのかしないのか？の考察が可能になる。 1)誤差を評価せよ。 1-A)コイルのインピーダンスの大きさにおける誤差の評価をせよ。 1-B)コイルの位相計算に含まれる誤差の評価をせよ。 1-C)コンデンサのインピーダンスの大きさにおける誤差の評価をせよ。 1-D)コンデンサの位相計算に含まれる誤差の評価をせよ。

課題解答例 1-A) コイルのインピーダンスの大きさにおける誤差の評価

[目的]

実験データに含まれる誤差について評価する。

[手順と結果]

変数の初期化する。

(%i1)

kill(all);

インピーダンスは以下の式で計算できる。

Impedance : Z=VA/VB *R

誤差は、以下で計算できる。

The Error : f=f(x0,x1,x2, ... xn)

df^2=(df/dx0*delt(x0))^2+(df/dx1*delt(x1))^2 ....+(df/dxn*delt(xn))^2

インピーダンス：Z を定義する。

(%i1)

Z:VA/VB*R;

各変数 VB,VA,R の微分項を計算する。

(%i2)

dVA:diff(Z,VA);

(%i3)

dVB:diff(Z,VB);

(%i4)

dR:diff(Z,R);

元の式で割ることで、相対誤差の計算を考える。

(%i5)

d1:dVA/Z;

(%i6)

d2:dVB/Z;

(%i7)

d3:dR/Z;

それぞれの誤差を、deltaVA,deltaVB,deltaR とする。

相対誤差は、次の式で計算できる。

(deltaZ/Z)^2 = (deltaVA/VA)^2 + (deltaVB/VB)^2+(deltaR/R)^2

以下の様に、相対誤差が定義できる。

(%i8)

deltaZZ:(d1*deltaVA)^2+(d2*deltaVB)^2+(d3*deltaR)^2;

式の簡単化を実行する。

(%i9)

expand(deltaZZ);

それぞれの変数について、相対誤差を仮定する。

VA は、1% 相対誤差として, deltaVA/VA=0.01、VB は、 1% 相対誤差として、 deltaVB/VB=0.01

最後に、R は、 10% 相対誤差として, deltaR/R=0.1 を代入する。

(%i10)

dZZ:ev(deltaZZ,deltaVB:0.01*VB,deltaVA:0.01*VA,deltaR:0.1*R);

相対誤差を％で求める。

(%i11)

dZZsq:sqrt(dZZ)*100;

測定データを入力する。周波数を datax に代入する。

(%i12)

datax:[500, 1000, 2000, 5000, 10000, 20000, 50000, 100000];

計算したインピーダンスを datay に代入する。

(%i13)

datay :[ 3.28, 6.20, 12.04, 30.22, 60.12, 121.63, 305.34, 709.08];

各値の誤差を計算する。

(%i14)

dy:datay*sqrt(dZZ);

別々のリストを [x1,y2], … [xn,yn] データ形式に変換する。

Ref.

[1] E. L. Woollett,"Maxima by Example Ch.1, Getting Started"

(%i15)

dataxy:map("[",datax,datay);

x-y データをグラフ表示する。

Ref.

[2] http://riotorto.users.sourceforge.net/gnuplot/errors/index.html

(%i16)

wxplot2d([discrete,dataxy]);

誤差を誤差棒で表示するために、package の draw を読み込む。

(%i17)

load(draw);

Fig.1 実験値のグラフ表示

x-y データリスト形式に変換する。

(%i18)

errdatay:map("[",datax,datay,dy);

実験値と誤差棒を同時に表示する。

(%i19)

wxdraw2d(

xrange=[200,200*10^3],

yrange=[1,1*10^3],

xlabel="f[Hz]",

ylabel="|Z|ohm",

logx=true,

logy=true,

grid =true,

error_type =y,

errors(errdatay),

/* 2nd graph for point*/

color = red,

point_size = 2,

point_type = circle,

points_joined = true,

points(dataxy)

);

Fig.2 実験値と誤差棒のグラフ表示

課題解答例 1-B) コイルの位相計算に含まれる誤差の評価について

[目的]

位相に含まれる誤差を評価して、グラフに表示する。

[手順と結果]

変数の初期化する。

(%i1)

kill(all);

位相の計算式を示す。

Phase : phi=360*dt*f

多変数における誤差の計算を以下に示す。

The Error : f=f(x0,x1,x2, ... xn)

df^2=(df/dx0*delt(x0))^2+(df/dx1*delt(x1))^2 ....+(df/dxn*delt(xn))^2

式を定義する。

(%i1)

phi:360*dt*f;

各変数の微分係数を計算する。

(%i2)

ddt:diff(phi,dt);

(%i3)

df:diff(phi,f);

変数の誤差を元の式で割り、相対値で表す。

(%i4)

d1:ddt/phi;

(%i5)

d2:df/phi;

各変数の誤差を deltadt, deltaf とする。

相対誤差は、以下で定義できる。

(deltaphi/phi)^2 = (deltadt/dt)^2 + (deltaf/f)^2

以上を定義する。

(%i6)

deltaphi:(d1*deltadt)^2+(d2*deltaf)^2;

式の簡単化を実行する。

(%i7)

expand(deltaphi);

値を代入する。時間 dt を 10% error にして、 deltadt/dt=0.1、

周波数 f が 1% error にして、 deltaf/f=0.01 とする。

(%i8)

dPhi:ev(deltaphi,deltadt:0.1*dt,deltaf:0.01*f);

相対誤差を%表示で求める。

(%i9)

dPhisq:sqrt(dPhi)*100;

次に、測定データを定義する。はじめに、周波数を datax に代入する。

(%i10)

datax:[500, 1000, 2000, 5000, 10000, 20000, 50000, 100000];

計算値の位相を datay に代入する。

(%i11)

datay :[ 75.6, 82.8, 86.4,90.0, 86.4, 89.3, 93.6, 90.0];

位相の誤差を計算する。

(%i12)

dy:datay*sqrt(dPhi);

実験データを[x1,y1] ……. [xn,yn]形式に変換する。

Ref.[1] E. L. Woollett,"Maxima by Example Ch.1, Getting Started"

(%i13)

dataxy:map("[",datax,datay);

実験データをグラフ表示する。

Ref.

[2] http://riotorto.users.sourceforge.net/gnuplot/errors/index.html

(%i14)

wxplot2d([discrete,dataxy]);

誤差棒を表示するために、package の draw を読み込む。

(%i15)

load(draw);

計算した誤差を、x-y 形式に変換する。

(%i16)

errdatay:map("[",datax,datay,dy);

Fig.1 実験データのグラフ表示

実験測定値と誤差棒を同時する。

(%i17)

wxdraw2d(

xrange=[100,200*10^3],

yrange=[50,100],

xlabel="f[Hz]",

ylabel="|Z|ohm",

logx=true,

grid =true,

error_type =y,

errors(errdatay),

/* 2nd graph for point*/

color = red,

point_size = 2,

point_type = circle,

points_joined = true,

points(dataxy)

);

Fig.2 実験データと誤差棒のグラフ表示

P

ha

se[

deg.

]