sanguisorbaのアイコン画像
sanguisorba 2021年05月07日作成 (2022年01月22日更新) © MIT
セットアップや使用方法 セットアップや使用方法 閲覧数 5960
sanguisorba 2021年05月07日作成 (2022年01月22日更新) © MIT セットアップや使用方法 セットアップや使用方法 閲覧数 5960

Wolfram Engine + Jupyter notebookで回路解析 (Windows)

Wolfram Mathematicaを個人で使うには高額すぎる。
そこで、無料で使えるWolfram EngineとJupyter notebookを使ってMathematica相当の環境を作ろうという話。

大学でWolfram Mathematicaのライセンスが無償配布されているという方はこの記事を見る必要ありません。

ちなみに、単純に回路解析をしたい場合はWolfram Alphaでも同じ事ができます。
https://www.wolframalpha.com/

後半の回路解析が間違っていたらコメントください

前書き

Q.どういうときにWolfram言語を使うの?
A.こういう回路を解析する時に使う。これを微積で計算するのは面倒くさい。

回路シミュレータによるとこのような出力電圧となるらしい。

1V 1000Hz

1V 7958Hz

Wolfram Engine 導入

大まかな流れ

  • ここから本体をダウンロード&インストール
    https://www.wolfram.com/engine/index.ja.php
  • Wolfram IDを作ってアクティベート(ダウンロード画面で次のステップとして書いてある事をするだけ)

これは大した作業じゃない。

Anaconda 導入

Jupyter notebookを用いるにはPythonの導入が必要ですが、Pythonを使った事がない人には敷居が高いです。
そこで今回はAnaconda Individual パッケージを使用します。このPythonパッケージにはJupyter notebookが予め内蔵されています。

プラットフォームに合わせてダウンロードしてください
https://www.anaconda.com/products/individual
これをインストーラの指示に従って入れてください。
・All Users でインストール
・Add Anaconda 3 to the system PATH environment variable のチェックは不要です。

プラグイン導入

このサイトから右上の緑色のCodeボタンを押してDownload Zipでダウンロードしてください。
https://github.com/WolframResearch/WolframLanguageForJupyter

解凍すると、WolframLanguageForJupyter-masterというフォルダの中にごちゃごちゃといろいろ入ってます。
このフォルダをそのままAnacondaのフォルダに移します。All Usersの場合、デフォルトのインストール先は
C:\ProgramData\Anaconda3
その後、Anaconda promptを起動して以下のコマンドを入力します

CD C:\ProgramData\Anaconda3\WolframLanguageForJupyter-master
.\configure-jupyter.wls add

その後、
jupyter kernelspec list
とコマンドを入力して、Available kernelsの中に"wolframlanguage"があれば成功です。
一回コマンドプロンプトを閉じて、改めてAnaconda 3 > Jupyter NotebookでJupyterを起動しましょう。

基本操作

Jupyter NotebookはWebベースのアプリケーションなので起動するとコマンドプロンプトと同時にブラウザが起動して、そこにJupyterの画面が表示されます。

右上のNewタブにWolfram Languageがあればちゃんと追加できています。これを選択しましょう。

キャプションを入力できます

すると次のような画面が表示されます。
キャプションを入力できます

In [] : のテキストボックスにWolfram言語を入力します。
例えば5+2を計算したい場合、5+2と入力して上のRunを押す。
例えば不定積分を行いたい場合は Integrate[Cos[x],x]と入力してRunを押す。
Wolram言語では三角関数の頭文字は大文字です。

そうすると次のように出力されます。
キャプションを入力できます

これがWolframです。

Wolfram言語については次のドキュメントを参照。
https://reference.wolfram.com/language/

前書きで書いた回路を分析する

コレ

今回は周波数特性をx軸に周波数、y軸に振幅比としてグラフに出力する事を目標とします。

Wolfram言語で回路解析をする場合はラプラス変換を用いて計算を行います。
そのため、回路方程式をラプラス変換によるs関数表記に変える必要があります。

回路を学習したことある人であれば、この3つの式は知っていると思います。

V=RIV=RI

V=1CIdtV=\frac{1}{C} \int I dt

V=LdIdtV=L\frac{dI}{dt}

s関数では、微分と積分を次のように表現します

1s=dt\frac{1}{s}= \int dt

s=ddts=\frac{d}{dt}

従って、コイルとコンデンサの電圧と電流の関係式を次のように書き換えます。

L[V]=1sCL[I]\mathcal{L}[V]=\frac{1}{sC}\mathcal{L}[I]

L[V]=sLL[I]\mathcal{L}[V]=sL\mathcal{L}[I]

これを基に回路方程式を立てます。

尚、ラプラス変換の記号は L\mathcal {L}、逆ラプラス変換の記号は L1\mathcal {L}^{-1}で表されます。

入力側の電流は

L[u(t)]=rL[I]+1sCL[I]=sCr+1sCL[I]\mathcal{L}[u(t)]=r\mathcal{L}[I]+\frac{1}{sC}\mathcal{L}[I]=\frac{sCr+1}{sC}\mathcal{L}[I]

L[I]=sCsCr+1L[u(t)]\mathcal{L}[I]=\frac{sC}{sCr+1} \mathcal{L}[u(t)]

出力側の電流は

L[y(t)]=RL[I]\mathcal{L}[y(t)]=R\mathcal{L}[I]

L[I]=1RL[y(t)]\mathcal{L}[I]=\frac{1}{R}\mathcal{L}[y(t)]

入力電流と反転出力電流が等しいと仮定して

sCsCr+1L[u(t)]+1RL[y(t)]=0\frac{sC}{sCr+1}\mathcal{L}[u(t)]+\frac{1}{R}\mathcal{L}[y(t)]=0

L[y(t)]=RsCsCr+1L[u(t)]\mathcal{L}[y(t)]=-R\frac{sC}{sCr+1}\mathcal{L}[u(t)]

一つ注意点として、入力電圧もラプラス変換をする必要があります。

ラプラス変換表から求めるのがセオリーですが、せっかくなのでWolframを用いて計算をしてみましょう。
LaplaceTransformの構文は

LaplaceTransform [**関数**, t, s]

です。中身の関数は時間tに関する表記となるように計算します。直流の場合は直流電圧をそのまま代入します。

以下は実行例です。(cos(ωt),5,eat,5cos(ωt)cos(\omega t), 5, e^{-at}, 5cos(\omega t))
キャプションを入力できます

前書きで載せたシミュレーションの入力電圧は交流1Vなので1×cos(ωt)1\times cos(\omega t)のラプラス変換した値である

L[u(t)]=ss2+ω2\mathcal{L}[u(t)]=\frac{s}{s^2+\omega ^2}

が入力電圧となります。ゆえに、

L[y(t)]=RsCsCr+1L[u(t)]=RsCsCr+1ss2+ω2\mathcal{L}[y(t)]=-R\frac{sC}{sCr+1}\mathcal{L}[u(t)]=-R\frac{sC}{sCr+1} \frac{s}{s^2+\omega ^2}

となります。
先程、ラプラス変換を行ったので今度はラプラス変換を行います。
さらに素子値を代入しても良いですが、この状態でラプラス変換を行う事もできます。
InverseLaplaceTransformの構文は

InverseLaplaceTransform [**関数**, s, t]

です。sとtが逆になります。具体的にはこんな感じ。

In [16] : InverseLaplaceTransform [-R*(s*C/(1+s*C*r))*(s/(s^2+ω^2)),s,t]

キャプションを入力できます

ちなみに、本家のWolfram Mathematicaで計算してもこんな感じ。
キャプションを入力できます

わかりやすく書くと、

y(t)=CRωsin(ωt)C2rRω2cos(ωt)C2r2ω2+1Rr(C2r2ω2+1)etCry(t)=\frac{CR\omega sin (\omega t) - C^2rR \omega ^2 cos (\omega t)}{C^2 r^2 \omega ^2 +1}-\frac{R}{r(C^2r^2 \omega ^2+1)} e^{- \frac{t}{Cr}}

これが入力電圧cosωtに対しての出力電圧です。ただし定常状態ではexp項がほぼ0となるので

y(t)=CRωsin(ωt)C2rRω2cos(ωt)C2r2ω2+1y(t)=\frac{CR\omega sin (\omega t) - C^2rR \omega ^2 cos (\omega t)}{C^2 r^2 \omega ^2 +1}

となります。

さっきの素子値を代入するとこんな感じになります。
r=200Ωr=200\Omega, R=100ΩR=100\Omega, C=100×109FC=100\times10^{-9}F

In [17] : InverseLaplaceTransform [-100*(s*10^(-7)/(1+s*10^(-7)*200))*(s/(s^2+ω^2)),s,t]

キャプションを入力できます
定常状態とすると、

y(t)=50000ωsin(ωt)ω2cos(ωt)5000000000+2ω2y(t)=\frac {50000\omega sin (\omega t)-\omega ^2 cos (\omega t)}{5 000 000 000 + 2 \omega ^2}

ただし、

ω=2πf\omega = 2 \pi f

これで振幅を求めてみます。NMaxValueで求めました。
キャプションを入力できます

1000Hzの時62.34mV、7958Hzの時353.56mVとなりました。
回路シミュレータでの値と比べてどうでしょう?
概ね一致しているのではないでしょうか。

カットオフ周波数は

f=12πfCrf=\frac{1}{2 \pi fCr}

です。シミュレーションで用いた素子値を代入すると7957.747となります。
この回路の最大振幅はRr=100200=0.5\frac{R}{r}=\frac{100}{200}=0.5です。
OPアンプでのカットオフ周波数は最大利得-3dBであり、振幅で表すと振幅Aに対してA2\frac{A}{\sqrt{2}}です。
よさげな計算結果が得られましたが、どうでしょう。
また、高周波になるにつれて最大振幅比の0.5,つまり500mVに近づきます。

In [34] : NMaxValue[(50000*(2*π*1000)*Sin[2*π*1000*t]-(2*π*1000)^2*Cos[2*π*1000*t])/(5000000000+2*(2*π*1000)^2),t]
Out[34] : 0.0623416

In [35] : NMaxValue[(50000*(2*π*7958)*Sin[2*π*7958*t]-(2*π*7958)^2*Cos[2*π*7958*t])/(5000000000+2*(2*π*7958)^2),t]
Out[35] : 0.353559

In [39] : NMaxValue[(50000*(2*π*7957.747)*Sin[2*π*7957.747*t]-(2*π*7957.747)^2*Cos[2*π*7957.747*t])/(5000000000+2*(2*π*7957.747)^2),t]
Out[39] : 0.353553

In [38] : 0.5/2^(0.5)
Out[38] : 0.353553

In [41] : NMaxValue[(50000*(2*π*50000)*Sin[2*π*50000*t]-(2*π*50000)^2*Cos[2*π*50000*t])/(5000000000+2*(2*π*50000)^2),t]
Out[41] : 0.493785

In [42] : NMaxValue[(50000*(2*π*100000)*Sin[2*π*100000*t]-(2*π*100000)^2*Cos[2*π*100000*t])/(5000000000+2*(2*π*100000)^2),t]
Out[42] : 0.498424

In [43] : NMaxValue[(50000*(2*π*200000)*Sin[2*π*200000*t]-(2*π*200000)^2*Cos[2*π*200000*t])/(5000000000+2*(2*π*200000)^2),t]
Out[43] : 0.499605

周波数と振幅比(最大値)をプロットすると、次のようになります。今回はPlotを用います。
キャプションを入力できます

入力はω-> 2πx, t -> yとして

In [51] : Plot [NMaxValue[(50000*2*π*x*Sin[2*π*x*y]-(2*π*x)^2*Cos[2*π*x*y])/(5000000000+2*(2*π*x)^2),y],{x,1,50000}, PlotRange ->{{0,50000},{0,0.5}}]

Plot内部でPlotRangeを指定してやらないとy軸の表示域が変な事になります。

前書きで書いた回路の出力電圧をラプラス変換無しで解く方法

コレ
今回提示した回路はコンピュータを使わないと解けないほど複雑という訳ではないです。ただ計算が面倒なだけ。

まず前提として覚えておかなければならない式

V=RIV=RI

V=1CIdtV=\frac{1}{C} \int I dt

V=LdIdtV=L\frac{dI}{dt}

これだけを用いて回路方程式を立てます。

まず入力電圧を決める必要があります。シミュレーションと同じ1Vの正弦波cos(ωt)cos (\omega t)とします。

入力側の電流は

rI+1CIdt=cos(ωt)rI+\frac{1}{C}\int I dt=cos(\omega t)

電荷に関する式に変換して微分方程式を解くと、(ここが面倒くさい)

dQdt+QCr=cos(ωt)r\frac{dQ}{dt}+\frac{Q}{Cr}=\frac{cos(\omega t)}{r}

今回は未定係数法で求めます。
通常係数はc1,c2...c_1, c_2 ...という感じに表しますが、コンデンサと被るのでbbで表記すると

Q=b1etCr+b2cos(ωt)+b3sin(ωt)Q=b_1e^{-\frac{t}{Cr}}+b_2 cos (\omega t)+b_3 sin (\omega t)

b2b_2, b3b_3は特解を代入すればよい。

d[b2cos(ωt)+b3sin(ωt)]dt+b2cos(ωt)+b3sin(ωt)Cr=cos(ωt)r\frac{d[b_2 cos (\omega t)+b_3 sin (\omega t)]}{dt}+\frac{b_2 cos (\omega t)+b_3 sin (\omega t)}{Cr}=\frac{cos(\omega t)}{r}

[b2ωsin(ωt)+b3ωcos(ωt)]+b2cos(ωt)+b3sin(ωt)Cr=cos(ωt)r[-b_2 \omega sin (\omega t)+b_3 \omega cos (\omega t)]+\frac{b_2 cos (\omega t)+b_3 sin (\omega t)}{Cr}=\frac{cos(\omega t)}{r}

ここからできる連立方程式

b2ω+b3Cr=0-b_2 \omega+\frac{b_3}{Cr}=0

b3ω+b2Cr=1rb_3 \omega+\frac{b_2}{Cr}=\frac{1}{r}

これを解いて

b2=Cω2C2r2+1b_2=\frac{C}{\omega ^2 C^2 r^2+1}

b3=ωC2rω2C2r2+1b_3=\frac{\omega C^2r}{\omega ^2 C^2 r^2+1}

だから、

Q=b1etCr+Cω2C2r2+1cos(ωt)+ωC2rω2C2r2+1sin(ωt)Q=b_1e^{-\frac{t}{Cr}}+\frac{C}{\omega ^2 C^2 r^2+1}cos (\omega t)+\frac{\omega C^2r}{\omega ^2 C^2 r^2+1}sin (\omega t)

b1b_1はt=0におけるコンデンサの電荷の初期値はQ=0だから

0=b1+Cω2C2r2+10=b_1+\frac{C}{\omega ^2 C^2 r^2+1}

b1=Cω2C2r2+1b_1=-\frac{C}{\omega ^2 C^2 r^2+1}

よって

Q=Cω2C2r2+1etCr+Cω2C2r2+1cos(ωt)+ωC2rω2C2r2+1sin(ωt)Q=-\frac{C}{\omega ^2 C^2 r^2+1}e^{-\frac{t}{Cr}}+\frac{C}{\omega ^2 C^2 r^2+1}cos (\omega t)+\frac{\omega C^2r}{\omega ^2 C^2 r^2+1}sin (\omega t)

あとは電荷を微分すれば電流になるので

I=dQdt=ω2C2rcos(ωt)ωCsin(ωt)ω2C2r2+1+1ω2C2r2+11retCrI=\frac{dQ}{dt}=\frac{\omega ^2 C^2 r cos (\omega t) - \omega C sin (\omega t)}{\omega^2 C^2 r^2 +1} + \frac{1}{\omega^2 C^2 r^2 +1} \frac{1}{r}e^{-\frac{t}{Cr}}

出力側の電流は

y(t)=RIy(t)=RI

I=y(t)RI=\frac{y(t)}{R}

入力電流と反転出力電流が等しいと仮定して

ω2C2rcos(ωt)ωCsin(ωt)ω2C2r2+1+1ω2C2r2+11retCr+y(t)R=0\frac{\omega ^2 C^2 r cos (\omega t) - \omega C sin (\omega t)}{\omega^2 C^2 r^2 +1} + \frac{1}{\omega^2 C^2 r^2 +1} \frac{1}{r}e^{-\frac{t}{Cr}}+\frac{y(t)}{R}=0

y(t)=ω2C2rRcos(ωt)+ωCRsin(ωt)ω2C2r2+11ω2C2r2+1RretCry(t)=\frac{-\omega ^2 C^2 r Rcos (\omega t) + \omega C Rsin (\omega t)}{\omega^2 C^2 r^2 +1} - \frac{1}{\omega^2 C^2 r^2 +1} \frac{R}{r}e^{-\frac{t}{Cr}}

これは先程の逆ラプラス変換で求めたものと一致しています。

y(t)=CRωsin(ωt)C2rRω2cos(ωt)C2r2ω2+1Rr(C2r2ω2+1)etCry(t)=\frac{CR\omega sin (\omega t) - C^2rR \omega ^2 cos (\omega t)}{C^2 r^2 \omega ^2 +1}-\frac{R}{r(C^2r^2 \omega ^2+1)} e^{- \frac{t}{Cr}}

ラプラス変換に頼らない回路計算がしたい方は、こちらをご参照ください。
技術者のための高等数学1 常微分方程式

sanguisorbaのアイコン画像
マイコンを使わない低レベルな電子工作とかPCBパターン製作など
ログインしてコメントを投稿する