0
$\begingroup$

My problem regards solving a differential equation and can be reduced to the problem of finding $f(x)$ such that $\frac{df}{dx}=\frac{dh}{df}$, with $h(f)$ a known function.

I have this list of data, from which I define a function myfun via Interpolation and enlarge its domain with Piecewise in the following way:

mydata={{0.0000255111, 3.48715}, {0.0000497289, 3.48715}, {0.0000760032, 
  3.68666}, {0.000101685, 4.13102}, {0.000127819, 
  4.8565}, {0.000160668, 5.89031}, {0.000199871, 
  7.03839}, {0.000251762, 8.11573}, {0.000510537, 
  10.1906}, {0.000759276, 10.6694}, {0.00131224, 
  10.8689}, {0.00205196, 10.8689}, {0.0039503, 10.8689}, {0.00586421, 
  10.8939}, {0.00936259, 10.9687}, {0.0121416, 11.0684}, {0.0157784, 
  11.1682}, {0.0214132, 11.7097}, {0.0269805, 12.465}, {0.0329414, 
  13.1633}, {0.0438442, 14.3354}, {0.0542594, 15.1783}, {0.0848459, 
  16.9373}, {0.105548, 18.4687}, {0.12379, 20.6365}, {0.132191, 
  22.4652}, {0.142875, 25.3083}, {0.15004, 27.8604}, {0.159849, 
  31.6752}, {0.167637, 34.0867}, {0.174909, 36.6097}, {0.180292, 
  39.116}, {0.187948, 41.5931}, {0.199271, 43.9}, {0.209739, 
  46.4811}, {0.225458, 48.9749}, {0.249185, 51.8844}, {0.270797, 
  54.2419}, {0.294283, 56.4441}, {0.320717, 58.4715}, {0.355922, 
  60.696}, {0.399908, 62.9903}, {0.468012, 65.264}, {0.537837, 
  67.6701}, {0.651061, 70.2674}, {0.770893, 72.5252}, {0.936818, 
  74.8075}, {1.12961, 76.9189}, {1.39071, 79.1248}, {1.74812, 
  80.9415}, {2.24823, 82.6618}, {2.80845, 83.7092}, {3.61942, 
  84.6744}, {4.4099, 85.2231}, {5.69251, 85.5709}, {7.33036, 
  85.8964}, {10.4519, 86.1783}, {13.1381, 86.691}, {16.5145, 
  87.616}, {20.7588, 88.7737}, {26.0938, 90.4665}, {31.4639, 
  92.2634}, {37.1584, 94.0748}, {44.8055, 96.3118}, {54.0264, 
  98.3997}, {66.5137, 100.666}, {83.6078, 102.479}, {105.095, 
  103.924}, {140.871, 105.274}, {176.744, 105.871}, {222.168, 
  106.231}, {279.265, 106.498}, {351.037, 106.642}, {441.253, 
  106.639}, {554.656, 106.842}, {697.203, 106.887}, {870.332, 
  106.902}, {1101.62, 106.932}, {1384.73, 106.932}, {1740.61, 
  106.932}, {2187.95, 106.932}, {2750.26, 106.932}, {3457.08, 
  106.932}, {4345.55, 106.932}, {5462.36, 106.932}, {6866.2, 
  106.932}, {8630.82, 106.932}, {10393.3, 106.932}, {12043.3, 
  106.932}};

myfun=Interpolation[mydata];
f[x_]:=Piecewise[{{myfun[x], x <= 1000}, {myfun[1000], x > 1000}}];

Now, I want to find a new function T[t] that satisfies this differential equation, where $f(T)$ is the function f defined with the list of data.

$$\dfrac{1}{T}\dfrac{dT}{dt}=-\dfrac{1}{1+\frac{1}{3}\frac{T}{f(T)}\frac{df(T)}{dT}}$$

This is what I've done (the initial conditions should be reasonable for the problem I am tackling):

diffeq={T[t]'/T[t]==-1/(1+0.3*T[t]/f[T[t]]*f[T[t]]')};
NDSolve=[diffeq,T[-40]==10^10,T,{t,-40,0}];

My code fails at finding the solution and I suspect that the problem resides in how I write the differential equation. I think that nesting the function T[t] that I want to find inside another function f[T] (and also making a derivative of f) drives Mathematica crazy.

How can I manage this task?

Thank you for your help!

Update, after the answer of Petrini

His method seems to work, but the error one receives is of the type

Power::infy: Infinite expression 1/0 encountered.

I have noticed that the function f has some serious smoothness problems, a glitch of Interpolation. Its shape is like this:

Graph of the function

but when you zoom in several intervals, like 10^-3,10^-2 you get this wiggling behaviour, whereas the data show a decreasing one, without oscillating.

function zoomed

The real problem (linked to the previous one I suppose) regards the derivative function, because if I plot f'[x] I get something that is a complete mess:

derivative of function

Hence, the 1/0 error should rise from here.

So maybe the problem now is more like

How do you get rid of that problematic behaviour of the function and its derivative?

Improve this question
$\endgroup$

2 Answers 2

Reset to default
2
$\begingroup$

This is an extended comment versus an answer, but you can improve your interpolation by performing it on log-log transformed. If you ListPlot mydata, you will see that much of the data is up against the x=0 axis.

ListPlot[mydata]

List Plot of mydata

Since your data looks much nicer and spread out on a ListLogLogPlot, you can perform the interpolation on Log-Log transformed data as shown:

lmyfun = Interpolation[Log@mydata];
lf[x_] := 
  Piecewise[{{Exp@lmyfun[Log[x]], x <= 1000}, {Exp@lmyfun[Log[1000]], 
     x > 1000}}];

Now, the interpolation is better behaved in the region you showed although you still have some discontinuities.

Show[{LogLogPlot[{f[x], lf[x]}, {x, 0.001, 0.01}, 
   PlotRange -> {{0.001, 0.01}, {10.6, 10.9}}, 
   PlotLegends -> {"original", "log transformed"}], 
  ListLogLogPlot[mydata]}]

Log transformed plot

Update to Include BSplineFunction

Interpolation will create a continuous function that passes through each interpolated point, but there is no guarantee that the derivatives will be continuous. As @Petrini stated in the comments, you may want to try fitting to have continuous derivatives. Fitting will not necessarily pass through data points, but it should have derivatives that are continuous up to a specified order. I will show an example using BSPlineFunction where I create a fit, then I resample with many points to create a smoother interpolation.

bsf = BSplineFunction[Log@mydata];
resampleddata = bsf[#] & /@ Subdivide[0, 1, 2000];
lbsmyfun = Interpolation[resampleddata];
lbsf[x_] := 
  Piecewise[{{Exp@lbsmyfun[Log[x]], 
     x <= 1000}, {Exp@lbsmyfun[Log[1000]], x > 1000}}];
Show[{LogLogPlot[{f[x], lf[x], lbsf[x]}, {x, 0.0001, 0.01}, 
   PlotRange -> {{0.001, 0.01}, {10.5, 10.9}}, 
   PlotLegends -> {"original", "log transformed", "bspline"}], 
  ListLogLogPlot[mydata]}]

BSpline Fit

Now, the derivatives are finite and positive so they show up on a log log plot without the breaks the OP observed.

LogLogPlot[lbsf'[x], {x, 0.0001, 1000}]

Derivative Plot

Update 2: Eliminate Piecewise Function by Padding Data

You have a division by zero warning due to the Piecewise function having a discontinuous derivative and extrapolation warnings due to mydata not going to low enough T values. You can eliminate these warnings by padding both ends by the initial and final values, respectively. You will need to make a decision whether it is valid to do so in this case.

mydata2 = 
  Table[{0.0000255111/2^(11 - i), 3.48715}, {i, 1, 10}]~Join~mydata~
   Join~Table[{12043.3 2^i, 106.932`}, {i, 1, 20}];
bsf = BSplineFunction[Log@mydata2];
resampleddata = bsf[#] & /@ Subdivide[0, 1, 4000];
lbsmyfun = Interpolation[resampleddata];
Plot[lbsmyfun[x], {x, -17, 23}]
lbsf[x_] := Exp@lbsmyfun[Log[x]]
diffeq = {T'[t]/T[t] == -1/(1 + (1/3)*T[t]/lbsf[T[t]]*lbsf'[T[t]])};
boundary = {T[-40] == 10^10};
odeqs = diffeq~Join~boundary;
tsol = NDSolveValue[odeqs, T, {t, -40, 0}];
LogPlot[tsol[t], {t, -40, 0}, PlotRange -> All]

Padded interpolation function

Improve this answer
$\endgroup$
7
  • $\begingroup$ This is useful, thank you. However, there still remain some discontinuities... It sounds really crazy that Mathematica is not able to give a simple interpolation, without those sharp oscillations! $\endgroup$
    – Lele
    Commented Oct 19, 2019 at 15:41
  • $\begingroup$ I added a BSplineFunction fit example. As suggested by @Petrini, you can have continuous derivatives with a fit. Interpolation goes through every point, but the derivatives may not be continuous. $\endgroup$
    – Tim Laska
    Commented Oct 19, 2019 at 19:16
  • $\begingroup$ Oh, very nice update! However, if I try to solve the differential equation I still obtain the 1/0 error. Damn... where does it come from? If I try with a domain like {-40,-25} it does not complain. But already with {-40,-20} I get the error... $\endgroup$
    – Lele
    Commented Oct 19, 2019 at 20:08
  • 1
    $\begingroup$ If you perform a FindRoot where T=1000, then you will see it occurs at t=-23.8819. Your Piecewise function does not have a continuous derivative @T=1000. You probably could pad your data out to 10^10 and eliminate the Piecewise function. $\endgroup$
    – Tim Laska
    Commented Oct 19, 2019 at 20:43
  • $\begingroup$ Ok, are you saying that prolonging the data imposing the continuity of the function in T=1000 makes it sharp? Do you have any other idea about how to smoothly pad the data? Than you! $\endgroup$
    – Lele
    Commented Oct 19, 2019 at 21:09
1
$\begingroup$

I could get some results by defining before the my data block:

diffeq = {T'[t]/T[t] == -1/(1 + (1/3)*T[t]/f[T[t]]*f'[T[t]])}
boundary = {T[-40] == 10^10}
odeqs = diffeq~Join~boundary

Where I took freedom to change 0.3 to (1/3), since they're not equal and your equation seems to use the latter.

And then using after defining f[x_]:

sol = NDSolve[odeqs, T, {t, -40, 0}];
Plot[Evaluate[{T[t] /. sol}], {t, -40, 0}]

But it found 1/0 values and supressed some of the evaluation. I hope it is of any help.

Improve this answer
$\endgroup$
3
  • $\begingroup$ Interesting! I get the 1/0 error too and maybe I know where the problem is, I post an update of the question! $\endgroup$
    – Lele
    Commented Oct 18, 2019 at 7:58
  • $\begingroup$ Maybe you could use fitting, instead of interpolation of the data. Therefore you can guarantee the smoothness of the derivative $\endgroup$
    – Petrini
    Commented Oct 18, 2019 at 12:49
  • $\begingroup$ Do you have any theoretical model for the data? You can fit into that $\endgroup$
    – Petrini
    Commented Oct 18, 2019 at 12:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.