# Fit a three-pulse sequence

I am wondering how to properly fit the three-pulse sequence shown below.

I have tried a method which is fitting the three pulses with three polynominal functions seperately and then finding the x-axis values where y==0. Afterwards, I can construct a piecewise function which contains the three polynominal functions with the boundary values same as the obtained x-axis values where y==0. The regime out of the boundary values is set to be 0. The fitting actually looks not too bad, as shown below.

However, as you may see, due to the inaccuracy of the set boundary values (the obtained x-axis values at y==0 are always approximate values), there are gaps such as at x=0.73 or x=1.5. Since I need to use this three-pulse sequence as a known function for NDSolve in the next stage, the gaps seem can cause errors like singularity or stiff system.

May I ask is there any better way to fit the three-pulse sequence continuously? The dataset can be found in the attached link.dataset

The fitting procedure what I have used is shown below:

PickedPumpData = Take[PumpData, {1, 49}];(*Data Pickup*)
RealPumpCurve = ListPlot[PickedPumpData, PlotRange -> All]
CurveFit = 6;
PumpFit[t_] = Fit[PickedPumpData, Table[t^i, {i, 1, CurveFit}], t]
SimulatedPumpCurve = Plot[PumpFit[t], {t, 0, 0.23}]
Solve[PumpFit[t] == 0, t]
Pumpcurve1[t_] := If[t > 0.2275646211603083 || t < 0, 0, PumpFit[t]]

PickedPumpData2 = Take[PumpData, {148, 191}];(*Data Pickup*)
RealPumpCurve2 = ListPlot[PickedPumpData2, PlotRange -> All]
CurveFit = 6;
PumpFit2[t_] = Fit[PickedPumpData2, Table[t^i, {i, 1, CurveFit}], t]
SimulatedPumpCurve2 = Plot[PumpFit2[t], {t, 0.735, 0.95}]
Solve[PumpFit2[t] == 0, t]
Pumpcurve2[t_] :=
If[t > 0.9477938759997302 || t < 0.73786902921016, 0, PumpFit2[t]]

PickedPumpData3 = Take[PumpData, {300, 340}];(*Data Pickup*)
RealPumpCurve3 = ListPlot[PickedPumpData3, PlotRange -> All]
CurveFit = 6;
PumpFit3[t_] = Fit[PickedPumpData3, Table[t^i, {i, 1, CurveFit}], t]
SimulatedPumpCurve3 = Plot[PumpFit3[t], {t, 1.5, 1.7}]
Solve[PumpFit3[t] == 0, t]
Pumpcurve3[t_] :=
If[t > 1.6922549659684891 || t < 1.4993056416906838, 0,
PumpFit3[t]]

WholePumpCurve[t_] := Pumpcurve1[t] + Pumpcurve2[t] + Pumpcurve3[t]
SimuWholePumpCurve = Plot[WholePumpCurve[t], {t, 0, 2}]

The following code is how the fitted function used in the NDSolve:

NrmlPumpFit[t] = A WholePumpCurve[t];
NSolve[Integrate[NrmlPumpFit[t], {t, 0, 2}] == PulseEnergy, A]
A = 18.994202212484236;
Nrmlpower = NrmlPumpFit[t]
Nrmlphoton = NrmlPumpFit[t]/(PlanckConst LightSpeed/(Wavelength))

OpticCoupl = 0.01;
optcrosection = 4 10^-17 (*cm^2*);
CrystalSize = 0.04189;
d = 0.3;
Ktot = 1/(9 10^-9)/10^3;
Kisc = 0.625 Ktot /10^3;
refrindex = 1.61;
Blight  = 3.87 10^21/10^3;;
ffactor = 0.5;
beamarea = Pi (1 10^-3)^2;
Wstimulatedemission = Nrmlpower 1.16 10;
w = 2 Pi 1.45 10^9;
q0 = 4277;
qthermal = 4277;
B = 9 10^-8 /10^3;
Kx = 2.7 10^4 /10^3;
Ky = 0.6 10^4 /10^3;
Kz = 1.7 10^3 /10^3;
Wxz = 1.1 10^4 /10^3;
Wyz = 2.2 10^4 /10^3;
Wxy = 0.44 10^4 /10^3;
Px = 0.76;
Py = 0.16;
Pz = 0.08;
Q0 = 7 10^3;
Qe = 7 10^3;
reciprocalQmprefactor = 7.26351*10^-19;

sol = NDSolve[{S0'[
t] == -(Wstimulatedemission + Nrmlphoton) OpticCoupl (1 -
Exp[-S0[t] optcrosection d/CrystalSize]) + (Krad +
Wstimulatedemission OpticCoupl) S1[t] + Kx X[t] + Ky Y[t] +
Kz Z[t],
S1'[t] == (Wstimulatedemission + Nrmlphoton) OpticCoupl (1 -
Exp[-S0[t] optcrosection d/CrystalSize]) - (Krad +
Wstimulatedemission OpticCoupl) S1[t] - Kisc S1[t],
X'[t] ==
Kisc S1[t]  Px - X[t] Kx - (X[t] - Z[t]) Wxz - (X[t] - Y[t]) Wxy -
B (X[t] - Z[t]) q[t],
Y'[t] ==
Kisc S1[t] Py - Y[t] Ky - (Y[t] - Z[t]) Wyz - (Y[t] - X[t]) Wxy,
Z'[t] ==
Kisc S1[t] Pz - Z[t] Kz + (X[t] - Z[t]) Wxz + (Y[t] - Z[t]) Wyz +
B (X[t] - Z[t]) q[t],
q'[t] == -(w (1/Q0 + 1/Qe -
reciprocalQmprefactor (X[t] - Z[t]))) 10^-3 (q[t] -
qthermal) + B (X[t] - Z[t]) q[t], q[0] == q0,
X[0] == Y[0] == Z[0] == S1[0] == 0, S0[0] == Ns}, {X, Y, Z, q, S0,
S1}, {t, 0, 3}]

Thank you!

• Why not construct an interpolation function? Also, it would be helpful if you showed the Mathematica code you used.
– JimB
Jan 10, 2020 at 18:07
• @JimB, Hi, thanks a lot for your reply! I have added the Mathematica code. May I ask can the interpolation function used in the NDSolve part? Jan 10, 2020 at 19:16
• I don't know if using Interpolation will work with NDSolve. However your dataset on DropBox has 902 data points and from the recently added code only the first 340 points are to be used. But in those 340 data points there are 68 times with the same first coordinate. Interpolation only works with unique first coordinates.
– JimB
Jan 10, 2020 at 20:19
• @JimB Thanks for pointing out the duplication of the coordinates. I will correct it and try to see if I can use Interpolation for NDSolve. Jan 13, 2020 at 11:13