Without even looking at your data, just your question, you need to understand how any such optimization tool works. Pretend you are a blind person, tasked with finding the location of lowest altitude on the surface of the earth. Yes, I know, it happens to be deep under water in the Pacific ocean. I'll give you scuba gear.
ll you are allowed to do is to use a cane, to infer the local clope of the surface. Then you will walk downhill. I told you to use the scuba gear! Things will get wet.
But now, suppose you were started out in some unfortunate place. Perhaps the vicinity of the Dead Sea? Thus a local depression where as you walk downhill, while you may get wet, you will also get stuck in the wrong place?
The point is, you need to use intelligent starting values for any such solver. Start it in a bad (wildly wrong) place, and you should expect to get complete crapola as a result.
DISCLAIMER: No blind individuals werre harmed in this thought experiment.
Having said all that, now we can look at your data.
Sigh. Can I suggest you learn to work in better units? Numbers with hugely varying scales will always cause you numerical problems on a computer. So instead of working with numbers that are on the order of 10^-11 seconds, use picoseconds or at least nanoseconds. Since a picosecond is 10^-12 seconds, that seems perfect here.
If nanoseconds would make you happier, then use them instead.
Next, consider shifting/translating the model. Since your data starts at roughly 18, then translate the model to reflect that. Otherwise, you will be using exponentials that start at 18. That causes you to have strange values appear in the coefficient a.
Sorry though. I won't use the app here, as I cannot as easily show and explain the results.
ft = fittype('a*cos(2*pi*f*(x-x0)).*exp(-(x-x0)/tp)','indep','x','problem','x0')
ft
ft =
General model:
ft(a,f,tp,x0,x) = a*cos(2*pi*f*(x-x0)).*exp(-(x-x0)/tp)
Here, when I plot the data, I can infer intelligent start points for the parameters.
a is the value at x == x0. That should be on the order of 1.
f and tp are rate parameters. tp is an exponential decay. 10 should be reasonable.
The most sensitive one is f, which as you have written it, should be roughly the inverse of the period, so a frequency. We can think of f as the number of oscillations per picosecond. Since the period of oscillation seems to be roughly 2, 0.5 is an intelligent start point. To get a decent estimate of the period of oscillation, I'll use the intersections tool written by Doug Schwarz, as found on the file exchange.
xpint = intersections(xp,y,[0 40],[0 0])
xpint =
18.6739787975366
19.7766155889863
20.8497137648191
21.9085954944223
22.9711077650663
24.0312152815146
25.0899375680738
26.1465450240289
27.2047855420838
28.2608584039546
29.3180979049269
30.3759822488665
31.4383534867009
32.4978593885987
33.5605930233276
34.622453287097
35.6820730748939
36.7440919564783
37.8046527675794
diff(xpint(1:2:end))
ans =
2.17573496728251
2.12139400024718
2.11882980300753
2.11484797400998
2.11331236284309
2.12025558177397
2.12223953662676
2.12148005156622
2.12257969268551
The period seems to be vaguely constant, though perhaps the first period is off from the remainder.
1/median(diff(xpint(1:2:end)))
ans =
0.471388153206563
Time to fit a model.
mdl = fit(xp,y,ft,'start',[1,.5,10],'problem',18)
mdl =
General model:
mdl(x) = a*cos(2*pi*f*(x-x0))*exp(-(x-x0)/tp)
Coefficients (with 95
a = 0.9925 (0.9156, 1.07)
f = 0.4445 (0.4405, 0.4485)
tp = 3.367 (2.971, 3.764)
Problem parameters:
x0 = 18
That seems to fit. At least, when I plot it, things seem to be somewhat reaonable.
plot(mdl)
hold on
plot(xp,y)
Now getting reasonable. This model suffers from a flaw, in that it assumes a fixed period of oscillation in that decay. Instead, your oscillations seem to be also slowing down. You could repair that, using a time varying coefficient f. Perhaps f(xp) might take the form f0 - f1*xp.
ft = fittype('a*cos(2*pi*(f0 - f1*(x - x0)).*(x-x0)).*exp(-(x-x0)/tp)','indep','x','problem','x0')
ft =
General model:
ft(a,f0,f1,tp,x0,x) = a*cos(2*pi*(f0 - f1*(x - x0)).*(x-x0)).*exp(-(x-x0)/tp)
mdl = fit(xp,y,ft,'start',[1,.5,.01,10],'problem',18)
mdl =
General model:
mdl(x) = a*cos(2*pi*(f0 - f1*(x - x0)).*(x-x0)).*exp(-(x-x0)/tp)
Coefficients (with 95
a = 0.9983 (0.9463, 1.05)
f0 = 0.4179 (0.4129, 0.4229)
f1 = -0.005473 (-0.006309, -0.004638)
tp = 3.494 (3.218, 3.771)
Problem parameters:
x0 = 18
figure
plot(mdl)
hold on
plot(xp,y)
That was a bit of a hack in the model. It does better, but we could surely do better yet with some effort. The problem is your fit will be jerked around by that first period. Since the absolute numbers are the largest there, the fit will be most impacted by errors in that first period. That would suggest a simple weighted fit might be as good, where the weight on the latter part of the series is increased.
