Chapter 2 2016.10.14 版 Maxima を用いた LC のインピーダンス測定について [ 目的 ] 電気通信大学 先進理工学科の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) コイルの位相計算に含まれる誤差の評価について
2-1 コイルにおけるインピーダンスの大きさと周波数の関係 [ 目的 ] 周波数依存性のグラフを表示する [ 結果 ] 変数をすべて消去する (%i1) kill(all); 直列抵抗 :r が存在するコイルのインピーダンスを z とする Fig.1 計算する回路図 (%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);
周波数の関数として インピーダンスを定義する (%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; Fig.1 計算する回路図 変数を正値と定義する (%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);
周波数の関数を定義 (%i6) define(azl(f),%o5); 理論値を表示する (%i7) wxplot2d(azl(f),[f,100,50*10^3], [logx],[logy], [xlabel,"frq"], [ylabel,"azl"], [gnuplot_preamble,"set grid"]); Fig.2 理論グラフ 実験データをリスト型で入力する (%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];
実験データを表示する (%i10) wxplot2d([discrete,datax,datay],[x,100,100*10^3], [logx],[logy],[style,points],[color,red],[gnuplot_preamble,"set grid"]); 理論値と実験値のグラフを表示する Fig.3 実験値グラフ (%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.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"]); Fig.5 LCR 値のグラフ 最後に 実験データと理論値と 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], [gnuplot_preamble,"set grid"]);
Fig.6 理論 LCR 値 実験値グラフ [Maxima ファイルの加工について ] 以上の結果は Maxima ver.5.23 を使用しました この version での日本語入力について Export で html 出力後に Word を用いて編集しています
2-3 コイルにおけるインピーダンスの位相と周波数の関係 [ 目的 ] 位相の周波数依存性グラフを表示する [ 結果 ] すべての変数を初期化する (%i1) kill(all); インピーダンス z を定義する (%i1) z:r+%i*w*l; Fig.1 計算する回路図 変数を正の値と定義する (%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 を周波数 :f に変換する (%i5) ratsubst(2*%pi*f,w,%o4);
周波数 f の関数として 位相を定義する (%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; Fig.1 計算する回路図 変数は 正の値とする (%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);
理論線を表示する (%i7) wxplot2d(thetal(f)*180/%pi,[f,100,50*10^3], [logx], [xlabel,"frq"], [ylabel,"azl"], [gnuplot_preamble,"set grid"]); 実験データをリストで定義して入力する Data from 1313068,sagara Fig.2 計算結果のグラフ (%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"]); 理論値と実験値を表示する Fig.3 実験値のグラフ (%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.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); 変数は 正の値と仮定する Fig.1 計算する回路図 (%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);
周波数の関数として インピーダンスを定義する (%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); Fig.1 計算する回路図 複素数の大きさを求める関数を使用 (%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);
周波数の関数を定義 (%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 理論グラフ 実験データをリスト型で入力する (%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];
実験データを表示する (%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 実験値グラフ Fig.4 理論と実験値グラフ
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); インピーダンス z を定義する (%i1) z:r+1/(%i*w*c); Fig.1 計算する回路図 変数を正の値と定義する (%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 を周波数 :f に変換する (%i5) ratsubst(2*%pi*f,w,%o4);
周波数 f の関数として 位相を定義する (%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); Fig.1 計算する回路図 複素数の大きさを求める関数を使用 (%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);
周波数の関数を定義 (%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"]); 実験データをリスト型で入力する Fig.2 理論グラフ (%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];
実験データを表示する (%i10) wxplot2d([discrete,datax,datay],[x,100,1*10^6], [logx],[style,points],[color,red],[gnuplot_preamble,"set grid"]); 理論値と実験値のグラフを表示する Fig.3 実験値グラフ (%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 理論と実験値グラフ
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]); Fig.1 実験値のグラフ表示 誤差を誤差棒で表示するために package の draw を読み込む (%i17) load(draw);
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]); Fig.1 実験データのグラフ表示 誤差棒を表示するために package の draw を読み込む (%i15) load(draw); 計算した誤差を x-y 形式に変換する (%i16) errdatay:map("[",datax,datay,dy);
Phase[deg.] 実験測定値と誤差棒を同時する (%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 実験データと誤差棒のグラフ表示