非線形フィッティング

測定データ

周波数:

> fdata:=map(`*`,[276.95, 281., 282.350, 282.950, 286.250, 286.9, 287.8, 288.7, 289.3, 291.150, 292.95],10^3);

fdata := [276950.00, 281000., 282350.000, 282950.00...
fdata := [276950.00, 281000., 282350.000, 282950.00...

インピーダンス:

> zdata:=[170.641, 57.0265, 15.3991, 45.5443, 1061.97, 4044.37, 633.092, 106.976, 337.904, 9660.67, 1387.93];

zdata := [170.641, 57.0265, 15.3991, 45.5443, 1061....
zdata := [170.641, 57.0265, 15.3991, 45.5443, 1061....

フィッティング:

> with(genfit);

[mnlfit]

> nda:=nops(fdata):
x := matrix(nda,1,fdata):
y := vector(nda,zdata):
w := vector(nda,[seq(1/j,j=1..nda)]):

> ini:=[1.17868*10^(-9),72.9823*10^(-12),4.35719*10^(-3),16.4963,72.9823*10^(-13),4.35719*10^(-2),16.4963];

ini := [.1178680000e-8, .7298230000e-10, .435719000...
ini := [.1178680000e-8, .7298230000e-10, .435719000...

> inits:=vector(7,ini);

inits := vector([.1178680000e-8, .7298230000e-10, ....
inits := vector([.1178680000e-8, .7298230000e-10, ....

> params:=[mnlfit(Zh,f,x,y,w,[Cb,Ca,La,Ra,Cc,Lc,Rc],inits,10^(-10))];

params := [.1536042361e-8, .7310282939e-10, .435543...
params := [.1536042361e-8, .7310282939e-10, .435543...

近似前と近似後を重ねてプロットします。

> para1:=convert(zip((a,b)->a=b,[Cb,Ca,La,Ra,Cc,Lc,Rc],params),set);

para1 := {Cb = .1536042361e-8, Ca = .7310282939e-10...
para1 := {Cb = .1536042361e-8, Ca = .7310282939e-10...

> new_Zh:=evalf(subs(para1, Zh)):

> para2:=convert(zip((a,b)->a=b,[Cb,Ca,La,Ra,Cc,Lc,Rc],ini),set);

para2 := {Cb = .1178680000e-8, Ca = .7298230000e-10...
para2 := {Cb = .1178680000e-8, Ca = .7298230000e-10...

> ini_Zh:=evalf(subs(para2, Zh)):

> data2:=zip((x,y)->[x,y], fdata, zdata);

data2 := [[276950.00, 170.641], [281000., 57.0265],...
data2 := [[276950.00, 170.641], [281000., 57.0265],...
data2 := [[276950.00, 170.641], [281000., 57.0265],...
data2 := [[276950.00, 170.641], [281000., 57.0265],...

プロットです。 ポイントは、測定データを表します。

赤い線がパラメータの初期値を使ったてプロットしたもの。青い線が求まったパラメータを使ってプロットしたものです。

> p1:=plots[pointplot](data2):
p2:=plot(ini_Zh,f=270000..295000, color=red):
p3:=plot(new_Zh,f=270000..295000, color=blue):
plots[display]({p1,p2,p3},view=-10000..20000,axes=boxed);

[Maple Plot]

求まったパラメータの値:

> para1;

{Cb = .1536042361e-8, Ca = .7310282939e-10, La = .4...
{Cb = .1536042361e-8, Ca = .7310282939e-10, La = .4...

初期値:

> para2;

{Cb = .1178680000e-8, Ca = .7298230000e-10, La = .4...
{Cb = .1178680000e-8, Ca = .7298230000e-10, La = .4...

>