Re: Complex Ratio Polynomial Fit
- To: mathgroup at smc.vnet.net
- Subject: [mg17291] Re: [mg17216] Complex Ratio Polynomial Fit
- From: Daniel Lichtblau <danl>
- Date: Fri, 30 Apr 1999 23:22:25 -0400
- References: <199904260520.BAA14395@smc.vnet.net.>
- Sender: owner-wri-mathgroup at wolfram.com
Colin L C Fu wrote:
>
> Hiya,
>
> I was told that the add-on packages in Mathematica are useful to solve
> more complicated
> and particular problems. Thanks!
>
> I tried to use Calculus'Pade' and Statistic'Nonlinear' to solve my
> initial problem. i am nearly there. The problems that I have when I
> used
> 1) Caluculus'Pade' is that the input has to be a function instead of
> data points.
>
> 2) Statistic'Nonlinear' is that I wouldnt be able to fit complex ratio
> polynomial in the form of (a1*x^n + a2*x^(n-1) + ... + I * (b1*x^m +
> ... ) )/(c1*x^n + c2*x^(n-1))
>
> Please advise
>
> Regards
> Colin
Here are a couple of possibilities. I do not recommend the first because
it is not numerically stable and also gives a denominator with complex
coefficients.
I work with the data below in both cases.
data = Table[{j,Exp[I*j]}, {j,1.,8.}];
<<Calculus`Pade`
func = InterpolatingPolynomial[data,x];
In[120]:= InputForm[res1 = Pade[func, {x,0,4,3}]]
Out[120]//InputForm=
(1.467035949378472 + 0.5407441618127256*I -
(0.859543402089528 + 0.35312868885585924*I)*x +
(0.6353138184138063 + 1.3596276434552435*I)*x^2 -
(0.5512404053403563 + 0.5305184226688553*I)*x^3 +
(0.10282711462013089 + 0.040427078567313424*I)*x^4)/
(1 + (0.27666936858971924 - 0.07097135892722134*I)*x +
(0.03937820386236435 - 0.023225588224089655*I)*x^2 +
(0.002734612910309003 - 0.003020589516318158*I)*x^3)
Second method. I did it by hand but one can use a Module, allow
different sizes for numerator and denominator in an obvious way.
avars = Array[a,5];
bvars = Array[b,5];
cvars = Array[c,4];
numer = (avars . x^(Range[Length[avars]]-1) +
I*bvars . x^(Range[Length[bvars]]-1));
denom = cvars . x^(Range[Length[cvars]]-1);
frac = numer/denom;
vars = Join[avars,bvars,cvars];
newdata = Map[(#[[2]]*denom - numer) /. x->#[[1]] &, data];
expr = Chop[Expand[newdata . ComplexExpand[Conjugate[newdata]]]];
{min, reps} = FindMinimum[Evaluate[expr],
Evaluate[Apply[Sequence, Map[{#,Random[]}&,vars]]],
MaxIterations->100];
In[122]:= InputForm[res2 = frac /. reps]
Out[122]//InputForm=
(0.02058985004730265 + 0.07929593673459871*x - 0.2062017221189517*x^2 +
0.06195974200201183*x^3 - 0.004726997214427548*x^4 +
I*(0.2170676356131192 + 0.5723762330590472*x - 0.2861909204202757*x^2
+
0.033511609426952804*x^3 - 0.0005167763595259002*x^4))/
(0.6344896295422635 + 0.22867115127738863*x - 0.1535436813404554*x^2 +
0.017373500484052767*x^3)
Below gives some idea of the quality of the results.
In[125]:= InputForm[Table[(res1 /. x->j) - data[[j,2]],
{j,Length[data]}]]
Out[125]//InputForm=
{0.00003703558987533384 - 0.000026321414304031343*I,
0.006330322722563231 - 0.003479356840526804*I,
0.10522966748468165 - 0.04231569229047634*I,
0.6636508460372128 - 0.17496932678184896*I,
2.434235674561071 - 0.3116331945135754*I,
6.275463430804625 + 0.06990806340173455*I,
12.567636872113988 + 1.9909110292933003*I,
20.807477504090674 + 6.4534166473977415*I}
In[126]:= InputForm[Table[(res2 /. x->j) - data[[j,2]],
{j,Length[data]}]]
Out[126]//InputForm=
{-0.6078178838576589 - 0.10384414156339561*I,
0.05033002569124129 - 0.13595018353257182*I,
0.2362133654530405 + 0.40201373487869096*I,
-0.3549795173154958 + 0.46306489695267666*I,
0.1906531107937608 - 0.9326975842276298*I,
0.4655235000902399 - 0.07799140240809915*I,
-0.19557262044268076 + 0.020032361403283128*I,
0.027803702327371435 + 0.0033266951709881587*I}
As a rule, your chances of getting a response will improve if you take
the trouble to type the package names correctly. Especially when you
mean "Statistics" rather than "Statistic".
Daniel Lichtblau
Wolfram Research