Message Boards Message Boards

GROUPS:

Obtain Prony Coefficients from the following data?

Posted 4 months ago
704 Views
|
11 Replies
|
1 Total Likes
|

I'm very new to Mathematica. I have been having the most difficult fitting my data to a Prony series to obtain coefficients. My data is listed as:

data ={{-7.41908,0.202581}, 
{-7.08513,0.22534}, 
{-6.80688,0.263415  },   
{-6.75203,0.27437   },    
{-6.47366,0.352769  },   
{-6.41908,0.372474  },   
{-6.14086,0.521425  },   
{   -6.13608,0.507203 },  
{   -6.08513,0.537803 },  
{   -5.80688,0.732571 },  
{   -5.8041,0.713165  },   
{   -5.75203,0.753371 },  
{   -5.47366,0.999272 },  
{   -5.4698,0.978897  },   
{   -5.41908,1.03047  },   
{   -5.14086,1.35962  },   
{   -5.13608,1.3339   },    
{   -5.08513,1.40366  },   
{   -4.80688,1.85884  },   
{   -4.8041,1.82376   },    
{   -4.75203,1.92504  },   
{   -4.47366,2.57341  },   
{   -4.4698,2.56778   },    
{   -4.41908,2.67567  },   
{   -4.36051,2.85621  },   
{   -4.14086,3.53946  },   
{   -4.13608,3.53409  },   
{   -4.08513,3.67948  },   
{   -4.02733,3.90283  },   
{   -3.80726,4.77778  },   
{   -3.80281,4.77227  },   
{   -3.75189,4.96328  },   
{   -3.69424,5.25015  },   
{   -3.47393,6.38345  },   
{   -3.46948,6.36833  },   
{   -3.41855,6.58298  },   
{   -3.36091,6.93619  },   
{   -3.1406,8.17887   },    
{   -3.13614,8.14628  },   
{   -3.02757,8.90897  },   
{   -2.80726,10.9395  },   
{   -2.80649,10.8011  },   
{   -2.80281,10.7765  },   
{   -2.69424,11.7196  },   
{   -2.47313,13.718   },    
{   -2.46948,13.6073  },   
{   -2.36091,14.8401  },   
{   -2.13981,17.9584  },   
{   -2.13614,17.8678  },   
{   -2.02757,19.8201  },   
{   -1.80648,25.1663  },   
{   -1.69424,28.8535  },   
{   -1.57108,34.0035  },   
{   -1.47315,39.0339  },   
{   -1.36091,46.0734  },   
{   -1.23773,55.7655  },   
{   -1.13981,65.2136  },   
{   -1.02757,78.5576  },   
{   -0.90441,96.9397  },   
{   -0.80648,114.854  },   
{   -0.69424,140.088  },   
{   -0.57108,174.772  },   
{   -0.47315,208.272  },   
{   -0.36091,254.819  },   
{   -0.23774,317.221  },   
{   -0.20832,328.147  },   
{   -0.13981,366.807  },   
{   0.095588,528.616  },   
{   0.125032,538.08   },    
{   0.19352,592.818   },    
{   0.428922,817.148  },   
{   0.458354,823.643  },   
{   0.526852,895.643  },   
{   0.762255, 1176.53 },  
{   0.791684, 1189.2  },   
{   0.860185, 1283.84 },  
{   1.095588, 1668.42 },  
{   1.125019, 1704.64 },  
{   1.19352,  1836.48  },   
{   1.221943, 1830.4  },   
{   1.226187, 1847.67 },  
{   1.428922, 1964.63 },  
{   1.458352, 1976.94 },  
{   1.555295, 2005.88 },  
{   1.55954,  2017.68  },   
{   1.762255, 2063.7  },   
{   1.791684, 2071.22 },  
{   1.888617, 2090.3  },   
{   1.892861, 2093.55 },  
{   2.095588, 2138.75 },  
{   2.125019, 2146.72 },  
{   2.221949, 2171.14 },  
{   2.226192, 2174.12 },  
{   2.428922, 2237.78 },  
{   2.458351, 2248.93 },  
{   2.555283  ,    2282.36    }, 
{   2.559527  ,    2285.98    }, 
{   2.791685  ,    2377.96    }, 
{   2.888616  ,    2407.15    }, 
{   2.892861  ,    2406.96    }, 
{   3.125019  ,    2492.99    }, 
{   3.221949  ,    2504.64    }, 
{   3.226192  ,    2504.83    }, 
{   3.458351  ,    2537.75    }, 
{   3.555282  ,    2527.52    }, 
{   3.559526  ,    2528.78    }, 
{   3.791685  ,    2494.84    }, 
{   3.888615  ,    2473.19    }, 
{   3.89286   , 2471.26 },  
{   4.221949  ,    2364.14    }, 
{   4.226192  ,    2359.41    }, 
{   4.555282  ,    2206.02    }, 
{   4.559526  ,    2193.85    }, 
{   4.888615  ,    1960.51    }, 
{   4.89286   , 1956.46 },  
{   5.221949  ,    1483.65    }, 
{   5.226192  ,    1474.91    }}

My prony series is listed as:

func=Normal [Series[((p_i^2)*E_i*(w^2))/(1+((p_i^2)*(w^2))),{w,0,22}]];

my fit is

FindFit[data,func,{E_i,p_i},w]

The results never fit the data.

I am very new to this and I am in great need of some assistance. I sure that the solution may be simple.

Thanks ahead of time.

11 Replies

Justin, your datapoints look like the overlay of two spectral lines. Why don't you try to fit them with two gaussians (given below as ffe)?

ffe = a Exp[-b (x - c)^2] + u Exp[-v (x - w)^2];
f2 = ffe /. FindFit[data,  ffe, {{a, 2000}, {b, 1}, {c, 4}, {u, 500}, {v, 1}, {w, 1.5}}, x]

and then

Plot[f2, {x, -2, 5}, Epilog -> Point /@ data]
Posted 4 months ago

Thank you Hans.

Fitting the results to a gaussian would be useful for the second half of these results. I currently need to fit them to an expansion that would get values to be eventually converted.

i borrowed your expression for Prony's equation. ffe = Ei+Normal[Series[(((pi^2)Ei(w^2))/(1+((pi^2)*(w^2)))),{w,0,19}]]

Where Prony =Ee+ Normal[Series[((((pi^2)Ei(w^2))/(1+((pi^2)*(w^2))))),{w,0,10}]] Pfit = Prony/. FindFit[data,Prony,{{Ee,1},{Ei,1},{pi,0.1}},w] SeriesCoefficent[Pfit,10]

gives me coefficents the results still dont fit. I may not be too sure how Prony's equation expands.

I am aware that people are able to do this in other software. Is there a way to convert this to sine, cosine?

The data plot is decreasing in both directions. Why do you expect a good fit to a Prony series?

Posted 4 months ago

My mistake :

data = {{   3.81E-08  ,    0.202581   }
{   8.22E-08  ,    0.22534    }
{   1.56E-07  ,    0.263415   }
{   1.77E-07  ,    0.27437    }
{   3.36E-07  ,    0.352769   }
{   3.81E-07  ,    0.372474   }
{   7.23E-07  ,    0.521425   }
{   7.31E-07  ,    0.507203   }
{   8.22E-07  ,    0.537803   }
{   1.56E-06  ,    0.732571   }
{   1.57E-06  ,    0.713165   }
{   1.77E-06  ,    0.753371   }
{   3.36E-06  ,    0.999272   }
{   3.39E-06  ,    0.978897   }
{   3.81E-06  ,    1.03047    }
{   7.23E-06  ,    1.35962    }
{   7.31E-06  ,    1.3339 }
{   8.22E-06  ,    1.40366    }
{   1.56E-05  ,    1.85884    }
{   1.57E-05  ,    1.82376    }
{   1.77E-05  ,    1.92504    }
{   3.36E-05  ,    2.57341    }
{   3.39E-05  ,    2.56778    }
{   3.81E-05  ,    2.67567    }
{   4.36E-05  ,    2.85621    }
{   7.23E-05  ,    3.53946    }
{   7.31E-05  ,    3.53409    }
{   8.22E-05  ,    3.67948    }
{   9.39E-05  ,    3.90283    }
{   0.000155861   , 4.77778 }
{   0.000157467   , 4.77227 }
{   0.000177057   , 4.96328 }
{   0.000202191   , 5.25015 }
{   0.000335792   , 6.38345 }
{   0.000339252   , 6.36833 }
{   0.000381459   , 6.58298 }
{   0.000435605   , 6.93619 }
{   0.000723441   , 8.17887 }
{   0.000730896   , 8.14628 }
{   0.000938485   , 8.90897 }
{   0.00155861    ,  10.9395  }
{   0.0015614 ,   10.8011   }
{   0.00157467    ,  10.7765  }
{   0.0020219 ,   11.7196   }
{   0.00336408    ,  13.718   }
{   0.00339252    ,  13.6073  }
{   0.00435605    ,  14.8401  }
{   0.00724748    ,  17.9584  }
{   0.00730896    ,  17.8678  }
{   0.00938484    ,  19.8201  }
{   0.0156141 ,   25.1663   }
{   0.020219  ,    28.8535    }
{   0.0268483 ,   34.0035   }
{   0.0336398 ,   39.0339   }
{   0.0435605 ,   46.0734   }
{   0.0578453 ,   55.7655   }
{   0.0724747 ,   65.2136   }
{   0.0938483 ,   78.5576   }
{   0.124621  ,    96.9397    }
{   0.156142  ,    114.854    }
{   0.20219   , 140.088 }
{   0.268485  ,    174.772    }
{   0.336397  ,    208.272    }
{   0.435605  ,    254.819    }
{   0.578436  ,    317.221    }
{   0.618985  ,    328.147    }
{   0.724746  ,    366.807    }
{   1.2462    ,  528.616  }
{   1.33362   , 538.08  }
{   1.56142   , 592.818 }
{   2.68486   , 817.148 }
{   2.87312   , 823.643 }
{   3.36397   , 895.643 }
{   5.78435   , 1176.53 }
{   6.18991   , 1189.2  }
{   7.24745   , 1283.84 }
{   12.462    ,  1668.42  }
{   13.3358   , 1704.64 }
{   15.6142   , 1836.48 }
{   16.6703   , 1830.4  }
{   16.834    ,  1847.67  }
{   26.8486   , 1964.63 }
{   28.7311   , 1976.94 }
{   35.9166   , 2005.88 }
{   36.2694   , 2017.68 }
{   57.8435   , 2063.7  }
{   61.8991   , 2071.22 }
{   77.3779   , 2090.3  }
{   78.1378   , 2093.55 }
{   124.62    ,  2138.75  }
{   133.358   , 2146.72 }
{   166.705   , 2171.14 }
{   168.342   , 2174.12 }
{   268.486   , 2237.78 }
{   287.31    ,  2248.93  }
{   359.156   , 2282.36 }
{   362.683   , 2285.98 }
{   618.992   , 2377.96 }
{   773.778   , 2407.15 }
{   781.377   , 2406.96 }
{   1333.58   , 2492.99 }
{   1667.05   , 2504.64 }
{   1683.42   , 2504.83 }
{   2873.1    ,  2537.75  }};

This master curve should not decrease on both sides. This was a mistake due to an aggregation of frequency curves.

I need to use Prony series to acquire the coefficents to plot the results in a time domain.

using E(t) = E0 +Ei*exp[-t/pi] (Summation of terms)

Put the data in proper Mathematica syntax. That means commas to separate list elements, and*10^(-8) instead of E-08.

Posted 4 months ago

Thanks Daniel:

its data = {{ 0.0000000381 , 0.2025810000 }, { 0.0000000822 , 0.2253400000 }, { 0.0000001560 , 0.2634150000 }, { 0.0000001770 , 0.2743700000 }, { 0.0000003360 , 0.3527690000 }, { 0.0000003810 , 0.3724740000 }, { 0.0000007230 , 0.5214250000 }, { 0.0000007310 , 0.5072030000 }, { 0.0000008220 , 0.5378030000 }, { 0.0000015600 , 0.7325710000 }, { 0.0000015700 , 0.7131650000 }, { 0.0000017700 , 0.7533710000 }, { 0.0000033600 , 0.9992720000 }, { 0.0000033900 , 0.9788970000 }, { 0.0000038100 , 1.0304700000 }, { 0.0000072300 , 1.3596200000 }, { 0.0000073100 , 1.3339000000 }, { 0.0000082200 , 1.4036600000 }, { 0.0000156000 , 1.8588400000 }, { 0.0000157000 , 1.8237600000 }, { 0.0000177000 , 1.9250400000 }, { 0.0000336000 , 2.5734100000 }, { 0.0000339000 , 2.5677800000 }, { 0.0000381000 , 2.6756700000 }, { 0.0000436000 , 2.8562100000 }, { 0.0000723000 , 3.5394600000 }, { 0.0000731000 , 3.5340900000 }, { 0.0000822000 , 3.6794800000 }, { 0.0000939000 , 3.9028300000 }, { 0.0001558610 , 4.7777800000 }, { 0.0001574670 , 4.7722700000 }, { 0.0001770570 , 4.9632800000 }, { 0.0002021910 , 5.2501500000 }, { 0.0003357920 , 6.3834500000 }, { 0.0003392520 , 6.3683300000 }, { 0.0003814590 , 6.5829800000 }, { 0.0004356050 , 6.9361900000 }, { 0.0007234410 , 8.1788700000 }, { 0.0007308960 , 8.1462800000 }, { 0.0009384850 , 8.9089700000 }, { 0.0015586100 , 10.9395000000 }, { 0.0015614000 , 10.8011000000 }, { 0.0015746700 , 10.7765000000 }, { 0.0020219000 , 11.7196000000 }, { 0.0033640800 , 13.7180000000 }, { 0.0033925200 , 13.6073000000 }, { 0.0043560500 , 14.8401000000 }, { 0.0072474800 , 17.9584000000 }, { 0.0073089600 , 17.8678000000 }, { 0.0093848400 , 19.8201000000 }, { 0.0156141000 , 25.1663000000 }, { 0.0202190000 , 28.8535000000 }, { 0.0268483000 , 34.0035000000 }, { 0.0336398000 , 39.0339000000 }, { 0.0435605000 , 46.0734000000 }, { 0.0578453000 , 55.7655000000 }, { 0.0724747000 , 65.2136000000 }, { 0.0938483000 , 78.5576000000 }, { 0.1246210000 , 96.9397000000 }, { 0.1561420000 , 114.8540000000 }, { 0.2021900000 , 140.0880000000 }, { 0.2684850000 , 174.7720000000 }, { 0.3363970000 , 208.2720000000 }, { 0.4356050000 , 254.8190000000 }, { 0.5784360000 , 317.2210000000 }, { 0.6189850000 , 328.1470000000 }, { 0.7247460000 , 366.8070000000 }, { 1.2462000000 , 528.6160000000 }, { 1.3336200000 , 538.0800000000 }, { 1.5614200000 , 592.8180000000 }, { 2.6848600000 , 817.1480000000 }, { 2.8731200000 , 823.6430000000 }, { 3.3639700000 , 895.6430000000 }, { 5.7843500000 , 1176.5300000000 }, { 6.1899100000 , 1189.2000000000 }, { 7.2474500000 , 1283.8400000000 }, { 12.4620000000 , 1668.4200000000 }, { 13.3358000000 , 1704.6400000000 }, { 15.6142000000 , 1836.4800000000 }, { 16.6703000000 , 1830.4000000000 }, { 16.8340000000 , 1847.6700000000 }, { 26.8486000000 , 1964.6300000000 }, { 28.7311000000 , 1976.9400000000 }, { 35.9166000000 , 2005.8800000000 }, { 36.2694000000 , 2017.6800000000 }, { 57.8435000000 , 2063.7000000000 }, { 61.8991000000 , 2071.2200000000 }, { 77.3779000000 , 2090.3000000000 }, { 78.1378000000 , 2093.5500000000 }, { 124.6200000000 , 2138.7500000000 }, { 133.3580000000 , 2146.7200000000 }, { 166.7050000000 , 2171.1400000000 }, { 168.3420000000 , 2174.1200000000 }, { 268.4860000000 , 2237.7800000000 }, { 287.3100000000 , 2248.9300000000 }, { 359.1560000000 , 2282.3600000000 }, { 362.6830000000 , 2285.9800000000 }, { 618.9920000000 , 2377.9600000000 }, { 773.7780000000 , 2407.1500000000 }, { 781.3770000000 , 2406.9600000000 }, { 1333.5800000000 , 2492.9900000000 }, { 1667.0500000000 , 2504.6400000000 }, { 1683.4200000000 , 2504.8300000000 }, { 2873.1000000000 , 2537.7500000000 }};

I'm still running into a bit of trouble. thanks for the continued help

Okay, this is usable data. I will make a few remarks and then show some possibly plausible code.

(1) It looks like the x axis values are almost logarithmic in scale. I threw them away and just worked with the y values. One can rescale the eventual result if need be.

(2) I still don't think a Prony series is the way to deliver a good fit.

(3) I do not have a good way to do this fit using only real valued exponent coefficients. I used Prony's method. It gives complex valued exponent coefficients. Code is cribbed from this MSE thread; there is similar somewhere on Wolfram Community as well.

(4) I also had to use trial and error to find a "good" cutoff for the number of terms. And I had to raise precision so that the numerics would behave.

(5) Even with all that, the fit is not good at the rightmost end, so I show the graph with the last few points removed. Imaginary parts also become nonzero, somewhat earlier. This is another sign that a Prony series might not be the best way to go.

vals = SetPrecision[data[[All, 2]], 200];
nn = Length[vals];

mat = Most[Partition[vals, Floor[nn/2], 1]];
svals = SingularValueList[mat, Tolerance -> 10^(-4)];
keep = Length[svals]
mat2 = Most[Partition[vals, keep, 1]];
rhs = Drop[vals, keep];
soln = PseudoInverse[mat2].rhs;

(* Out[162]= 30 *)

So our series will have 30 terms.

roots = xx /. NSolve[xx^keep - soln.xx^Range[0, keep - 1] == 0, xx];
freqs = Log[roots];
newmat = Map[roots^# &, Range[0, nn - 1]];
coeffs = PseudoInverse[newmat].vals;
newf = TrigExpand[ExpToTrig[coeffs.Exp[freqs*(t-1)]]];

(I offset the time parameter so that at t=1 we would get the value at zero.)

We now have a Prony series. We can plot how well (or poorly) it fits.

p1 = Plot[Re@newf, {t, 1, nn - 7}, PlotRange -> All, 
  WorkingPrecision -> Precision[newf]]

enter image description here

Now plot the data values, connected by segments.

p2 = ListLinePlot[data[[All, 2]], ColorFunction -> (Red &)]

enter image description here

Superposed:

Show[p1, p2]

enter image description here

It might be easier to see the fit if we do not connect the data points.

p2Points = ListPlot[Style[data[[All, 2]], Red]];
Show[p1, p2Points]

enter image description here

--- edit ---

This works better if the data is first smoothed. When I redo the computations I get no trouble at the right end, and no imaginary part, provided I use enough singular values (correspondiong to Prony series terms). I tried using both listConvolve and GaussianFilter. The latter worked slightly better so I'll show the results from that.

origvals = SetPrecision[data[[All, 2]], 200];
vals = GaussianFilter[origvals, 8, 
   WorkingPrecision -> Precision[origvals]];
nn = Length[vals];
mat = Most[Partition[vals, Floor[nn/2], 1]];
svals = SingularValueList[mat, Tolerance -> 10^(-5)];
keep = Length[svals]
mat2 = Most[Partition[vals, keep, 1]];
rhs = Drop[vals, keep];
soln = PseudoInverse[mat2].rhs;

roots = xx /. NSolve[xx^keep - soln.xx^Range[0, keep - 1] == 0, xx];
freqs = Log[roots];
newmat = Map[roots^# &, Range[0, nn - 1]];
coeffs = PseudoInverse[newmat].vals;
newf = TrigExpand[ExpToTrig[coeffs.Exp[freqs*(t - 1)]]];

p1 = Plot[newf, {t, 1, nn}, PlotRange -> All, 
   PlotStyle -> Directive[Thickness[.02], Dashed, Green], 
   WorkingPrecision -> Precision[newf]];
p2 = ListPlot[vals, PlotStyle -> Directive[Red, PointSize[.01]]];
Show[p1, p2]

enter image description here

--- end edit ---

Posted 3 months ago

Forgive me , but I'm still a little confused with results. How do you acquire the coefficients from this?

pi and Ei ?

What I computed is a bit different from your proposed form. It gives a sum of complex exponential. The coefficients are 'coeffs' and the exponents are 'freqs' in my code. Look up "Prony's method" for more details on what it provides.

Sorry for not knowing enough english: what is a Prony series? And what do you mean by

E(t) = E0 +Ei*exp[-t/pi] (Summation of terms) ?

What is Ei, and what is pi in your problem?

Perhaps you could write down the first 3 or four terms of your series, if possible in Mma-consistent form.

Concerning your new data given above you might want to have a look at

pp = FindFit[data, (a t)/(1 + b t), {a, b}, t]
f3 = (a t)/(1 + b t) /. pp
Plot[f3, {t, 0, 150}, Epilog -> Point@data, PlotRange -> All]

Ok. In the meantime I figured out that Prony seems to be a name.

Do you perhaps mean something like this (ignore the error message that appears - at least in my system)

prony = e0 + Sum[e[i] Exp[-(t/tau[i])], {i, 1, 19}]
param = Flatten[Join[{e0}, Table[{e[i], tau[i]}, {i, 1, 19}]]]
fit = FindFit[data, prony, param, t]
Plot[prony /. fit, {t, 0, 150}, Epilog -> Point /@ data,  PlotRange -> All]
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract