007_Reluctance.mw

非線形電気抵抗と微分方程式

イントロダクション

このワークシートは非線形の磁気抵抗をモデル化し、その解とグラフを求めます。

> restart;

微分方程式の定義

d1, d2, および d3を定義し、方程式を作成します。

> d1 := Vi(t) - x1(t)*r1 - l1*diff(x1(t),t) - N1 * diff(x0(t),t) = 0;

d1 := Vi(t)-x1(t)*r1-l1*(diff(x1(t), t))-N1*(diff(x0(t), t)) = 0

> d2 := 0 - x2(t)*(r2+R1) - (l2+L1)*diff(x2(t),t) - N2 * diff(x0(t),t) = 0;

d2 := -x2(t)*(r2+R1)-(l2+L1)*(diff(x2(t), t))-N2*(diff(x0(t), t)) = 0

> d3 := N1 * diff(x1(t),t) + N2 * diff(x2(t),t) = Rel * diff(x0(t),t);

d3 := N1*(diff(x1(t), t))+N2*(diff(x2(t), t)) = Rel*(diff(x0(t), t))

パラメータの定義

> Vi(t) := 650*cos(w*t):
w := 100*Pi:
N1 := 70:
N2 := 35:
r1 := .53:
r2 := .13:
l1 := .32*10^(-3):
l2 := .08*10^(-3):
R1 := 100:
L1 := 100*10^(-3):

非線形電気抵抗

v[x] および v[y] に対応するベクトルを定義します。

> vx := [-30*10^(-3),-25*10^(-3),-20*10^(-3),-15*10^(-3),-10*10^(-3),-5*10^(-3),0,5*10^(-3),10*10^(-3),15*10^(-3),20*10^(-3),25*10^(-3),30*10^(-3)];

vx := [(-3)/100, (-1)/40, (-1)/50, (-3)/200, (-1)/100, (-1)/200, 0, 1/200, 1/100, 3/200, 1/50, 1/40, 3/100]

> vy := [400000,160000,40000,23000,22000,21000,20000,21000,22000,23000,40000,160000,400000];

vy := [400000, 160000, 40000, 23000, 22000, 21000, 20000, 21000, 22000, 23000, 40000, 160000, 400000]

Maple のスプライン関数を利用し、スプラインを x[0] の関数として生成します。関数は区分分割関数として完全に解析的な形で表現されます。また、ここで Maple の Spline 関数では、関数 (ここでは x[0](t) ) を変数として扱うことができませんので、ひとまず x の関数と定義し、後で x を x[0](t) と置き換えます。

> Rel := spline(vx,vy,x,cubic);

Rel := PIECEWISE([815200000/193+659841600000/1351*x+24382080000000/1351*x^2+270912000000000/1351*x^3, x < (-1)/40], [573400000/1351+6268800000/193*x-256320000000/1351*x^2-57600000000000/1351*x^3, x < ...

> Rel := subs(x=x0(t), Rel);

Rel := PIECEWISE([815200000/193+659841600000/1351*x0(t)+24382080000000/1351*x0(t)^2+270912000000000/1351*x0(t)^3, x0(t) < (-1)/40], [573400000/1351+6268800000/193*x0(t)-256320000000/1351*x0(t)^2-57600...

解決方法

初期条件を追加します。

> inits := x0(0)=0, x1(0)=0, x2(0)=0;

inits := x0(0) = 0, x1(0) = 0, x2(0) = 0

問題を RK4 を使ってステップ・サイズ .0005で解きます。

> sol3 := dsolve({d1,d2,d3,inits},{x0(t),x1(t),x2(t)}, type=numeric);

sol3 := proc (x_rkf45) local res, data, vars, solnproc, outpoint, ndsol, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; _EnvDSNumericSaveDigits := Digits; Digits := 14; if...

に対応した解をプロットします。

> plots[odeplot](sol3,[t,x1(t)],0.. .05,numpoints=400);

[Plot]