0
$\begingroup$

We have two lists (vectors) of data, y and x, we can imagine x being time steps (0,1,2,...) and y some system property computed at each value of x. I'm interested in calculating the derivative of log of y with respect to log of x, and the question is how can I preform such calculations in Mathematica? We first go to log space, e.g. logx = Map[Log,x] similarly for y. Then what method do we use for the differentiation dlog(y)/dlog(x)? and then eventually we'll plot it as a function of log(x).

Context: When we have a quantity y for which we expect to see a power law behaviour as a function of x, so $y \propto x^{\alpha},$ taking the log derivative of y w.r.t to x allows us to extract the actual $\alpha$ corresponding to our data.

Improve this question
$\endgroup$
1
  • $\begingroup$ Could interpolate in log-log space and take derivatives of the interpolated form. This needs to be done with a reasonable interpolation degree so as not to introduce artificial wiggles. $\endgroup$
    – Daniel Lichtblau
    Feb 26, 2018 at 16:47

2 Answers 2

Reset to default
2
$\begingroup$

In principle, I like most of all the approach proposed by @Daniel Lichtblau as a comment. I would propose Daniel to formulate it as a regular answer.

If the data is noisy, the interpolation can give rise to problems, since it goes through all points. The function, thus, strongly oscillates, and its derivatives can appear unnecessarily large and also noisy.

Instead, the strategy using fitting may be a good idea to apply, if the noise is present. The latter enables one to create a smooth, slowly varying function, thus effectivey reducing the effect of the noise in the derivative. Here I would like to propose a simple idea of how to treat the problem by the fitting. Since the author gave no list, to be able to show code, I will create such a list. Since the author is primarily interested in a power law-type of dependence, this will be the cases of the lists below.

The first example is a simple power law to which I add a bit of noise:

lst1 = Table[{x, x^(3/2) + RandomReal[{-30, 30}]}, {x, 1, 100}];
ListPlot[lst1]

yielding this:

enter image description here

This list can be transformed into a Log-Log list:

lst1Log = lst1 /. {x_, y_} -> {Log[x], Log[y]} // N;

The next step is to fit it to a suitable function. In principle, by looking at its plot one can easily guess that at large values of t=Log[x] it has a linear asymptote. This agrees with the author's expectation of a power law. It is not astonishing, since I created this list to be such. However, the authors expectation and the form of the plot gives me a hint to fit it to a linear function:

     model1 = a + b*t;
ff1 = FindFit[Select[lst1Log, #[[1]] > 3 &], model1, {a, b}, t]
Show[{
  ListPlot[lst1Log, PlotStyle -> Blue],
  Plot[model1 /. ff1, {t, 0, 5}, PlotStyle -> Red]
  }]
    (*   {a -> 0.167, b -> 1.46}   *)

where a and b are the fitting parameters, bbeing the exponent of the power law. The result differs from the power law b=1.5 since it is influenced by stattered points at small t. To improve it one can cut away the list points at smaller t: Select[lst1Log, #[[1]] > 4 &] and make a fit starting from the value b=1.46 obtained within the previous step:

model1 = a + b*t; ff1 = FindFit[Select[lst1Log, #[1] > 4 &], model1, {a, {b, 1.46}}, t] Show[{ ListPlot[lst1Log, PlotStyle -> Blue], Plot[model1 /. ff1, {t, 0, 5}, PlotStyle -> Red] }]

(* {a -> -0.0245, b -> 1.51} *)

Here how it looks like:

enter image description here

The difference from 1.5 is due to the noise.

The second example is the list containing two power laws and exhibiting a slightly more complex dependence:

    lst2 = Table[{x, 5 x^(1/2) + 0.001*x^2.7 + RandomReal[{-10, 10}]}, {x,
         4, 1000}];
    lst2Log = lst2 /. {x_, y_} -> {Log[x], Log[y]} // N; 
ListPlot[lst2]

enter image description here

It is the same story, but looking at the plot of the Log-Log list (that is, lst2Log) one can guess that it can be fitted by a slightly different function:

model2 = (a + b*t^(d + 1))/(1 + c*t^d)

Here a, b, c and d are parameters. In such functions, however, the exponents (such as d) are usually difficult to fit. For this reason instead of the fitting it is convenient to manipulate with the parameter dwhile fitting all the rest:

Manipulate[
 model2 = (a + b*t^(d + 1))/(1 + c*t^d);
 ff2 = FindFit[lst2Log, model2, {{a, 3}, b, c}, t];
 Show[{
   ListPlot[lst2Log, PlotStyle -> {Blue, PointSize[0.008]}],
   Plot[model2 /. ff2, {t, 0, 7}, PlotStyle -> Red]
   }], 
{{d, 3.58}, 2, 5, Appearance -> "Labeled"}]

By manipulating I found that the value d=3.58 gives a good approximation. It is this value that I placed as an initial one for convenience. This is what you should see on the screen:

enter image description here

Now one can get the approximation more precisely, by using the values of the parameters obtained in the previous step as their initial values:

model2 = (a + b*t^(4.6 + eps))/(1 + c*t^(3.6 + eps));
ff2 = FindFit[Select[lst2Log, #[[1]] > 3.5 &], 
  model2, {{a, 2.55}, {b, 0.007}, {c, 0.003}, {eps, 0.2}}, t]
Show[{
  ListPlot[Select[lst2Log, #[[1]] > 3.5 &], 
   PlotStyle -> {Blue, PointSize[0.008]}],
  Plot[model2 /. ff2, {t, 0, 7}, PlotStyle -> Red]
  }]

(*  {a -> 2.72, b -> 0.00462, c -> 0.00228, eps -> 0.258}  *)

enter image description here

The exponent in infinity will in this case be b/c. However you can here find the exponent value in each point by calculating the derivative of the model2in the desired point.

Generally, in this approach the fitting gives you a simple smooth approximation of the desired function, the derivative of which you can further find.

Have fun!

Edit 1: To illustrate how to find the regimes. Within the example 2 you can define the derivative of ydenoted below as der[t] and plot it:

Clear[der];
der[t_] := Evaluate[D[model2 /. ff2, t] // Simplify];

Plot[der[t], {t, 0, 7}, PlotStyle -> Red, 
 AxesLabel -> {Style["t=Log[x]", Italic, 12], 
   Style["der=\!\(\*FractionBox[\(dLog[y]\), \(dt\)]\)", Italic, 12]}]

Here you will find 3 regimes, from 1 to several decades each. In the middle a regime with a derivative described by an almost straight tilted line is visible. It is flanked by two power law rerimes at small and large t values.

Edit 2: To address your question with your list.

In this case I would do as follows. Here is your list, that should be somehow downloaded

lst = Drop[
Import[FullPathToYourList, "Table"], 1] ;

I have dropped the first element of the list, since there is something strange in it. But it makes a negligible influence on the result. and lst1=lst/. {x_, y_} -> {Log[x], Log[y]} Now, looking at the ListPlot[lst1] one sees that it consists of two straight regions corresponding to two different regimes (i) at Log[x]<4 and (ii) at Log[x]>4. In principle this means that you could express it in the form of a formula y=a e^{-\frac{x}{d}} x^{\alpha }+b \left(1-e^{-\frac{x}{d}}\right) x^{\beta } (I hope that the TeX expression here works). Here d is a cut-off close to 10^4. However, it is not easy to fit such a formula and to do it on the huge interval that takes place here.

For this reason it is more easy to fit the log-log list lst1, and to do it for t=Log[x]<4´separately from that for t>4. Below I do it:

   model1 = a + b*x;
model2 = c + k*x;

ff1 = FindFit[Select[lst, #[[1]] < 4 &], model1, {a, b}, x]
ff2 = FindFit[Select[lst, #[[1]] > 10 &], model2, {c, k}, x]
Show[{
  ListPlot[lst, PlotStyle -> Blue, 
   AxesLabel -> {Style["t=Log[x]", 12, Italic], 
     Style["Log[y]", 12, Italic]}],
  Plot[model1 /. ff1, {x, 1, 4}, PlotStyle -> Red],
  Plot[model2 /. ff2, {x, 4, 14}, PlotStyle -> Green]
  }]

(*  {a -> 1.54653449015407020, b -> 0.760890314643084274}

{c -> 3.002027525233422865, k -> 0.3272780474976788419}   *)

On the screen you will see the following plot:

enter image description here

Now, the fitting parameters b and k above are the exponents you are looking for. So I would say that the exponents are 0.76 and 0.33.

Pay attention that in the ff2 I fit the part of the list with t>10. It is because the interval 4<t<10 looks as a transition part of the curve. It is for this reason the solid green line deviates from the blue points in this interval. This is, however, up to you to choose how and where it is better to fit. It strongly depends on the nature of the problem.

Have fun!

Improve this answer
$\endgroup$
14
  • $\begingroup$ @user929304 The answer is "yes". This is exactly the second example in my answer. If you can find there two different regions covering not less than a decade along t each, such that within each such a region you can accurately approximate the list by a straight line, you get different regimes of your system. This is easy to see by taking D[model2 /. ff2, t] and plotting it. Please have a look at the edit $\endgroup$
    – Alexei Boulbitch
    Feb 27, 2018 at 12:52
  • $\begingroup$ @user929304 The best would be, if you have a theory describing these different regimes. I think to have dealt with such a situation. You may have a look into papers Uhde et al. Phys. Rev. E 72, 061916 (2005); Uhde et al. Phys. Rev. Lett. 94, 198102 (2005); Ter-Oganessian et al. Phys. Rev. E 72, 041510 (2005) and Ter-Oganessian et al. Phys. Rev. E 72, 041511 (2005). There we also faced with a situations with three different power law regimes and the necessity to fit them all. But there the theory has been first created. $\endgroup$
    – Alexei Boulbitch
    Feb 27, 2018 at 13:03
  • $\begingroup$ I'm playing catch-up this week so probably won't have time to try code for this. Wanted to mention that a related possibility is to use a smoothing filter as preprocessor. It would be in conjunction with the various other tactics under consideration, not a replacement for them. $\endgroup$
    – Daniel Lichtblau
    Feb 27, 2018 at 17:10
  • $\begingroup$ @user929304 Please have a look at the last edit. $\endgroup$
    – Alexei Boulbitch
    Mar 5, 2018 at 9:41
  • $\begingroup$ @AlexeiBoulbitch Thank you very much, this has been incredibly instructive! $\endgroup$
    – user52181
    Mar 5, 2018 at 15:22
3
$\begingroup$

While it seems unlikely that this is what you will actually want, it is what you asked for:

xs = Range[100]
ys = 3*xs^0.5
Differences[Log@ys]/Differences[Log@xs]

Note that you cannot have a zero value in xs.

Improve this answer
$\endgroup$
2
  • $\begingroup$ Thanks a lot. This is definitely very close, I wonder, is there more accurate ways of performing the numerical differentiation? (for instance higher order differences). Additionally, when taking the Differences, the resulting list has one element less than the involved lists (xs e.g.,) so for plotting one element from xs has to be removed. $\endgroup$
    – user52181
    Feb 26, 2018 at 15:48
  • $\begingroup$ Under your specification, the relationship is loglinear, so in this case no, there is no more accurate method. In other circumstances, you might average forward and backwards differences or fit a parameterized function and use its derivatives. $\endgroup$
    – Alan
    Feb 26, 2018 at 17:56

Your Answer

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