3
$\begingroup$

I am using the Drubin numerical inversion out of Laplace space, and need to see how the function performs over a wide range of values, covering multiple logrithmic cycles. However, the inversion is slow, and when I need the definition at small times, it takes a very long time to run through the whole range. I tried the Crump accelerator, but it is very memory intensive, so I'd rather not use it.

The code I have at the moment is:

adurb[tmax_, alpha_, tol_] := adurb[tmax, alpha, tol] = alpha - Log[tol]/(8 tmax) // N
cdurb[F_, s_, tmax_, alpha_, tol_, k_] := F /. s -> adurb[tmax, alpha, tol] + I Pi k/(4 tmax) // N
cdurb[F_, s_, tmax_, alpha_, tol_, 0] := F/2 /. s -> adurb[tmax, alpha, tol] // N

NLInvDurb[F_, s_, t_, j_, tmax_, alpha_, tol_] := NLInvDurb[F, s, t, j, tmax, alpha, tol] = Exp[adurb[tmax, alpha, tol] t]/(4 tmax)*Sum[Re[cdurb[F, s, tmax, alpha, tol, k]] Cos[k Pi t/(4 tmax)] - Im[cdurb[F, s, tmax, alpha, tol, k]] Sin[k Pi t/(4 tmax)], {k, 0, j}] // N    

ini = 0.001;
imax = 1000;
istep = 0.01;
NStep = 1000;

datrul1 = {CD -> 0.1, CfD -> 16, ϕ -> 10, α -> 0.3, η -> 10, S -> 0.1};

qD1[u_] = 1/u^2 (1 + 
 u^2 CD (Pi/CfD 1/(u (u/η + (2 Sqrt[u])/CfD)^(1/2))) + 
 u S CD)/((1 + ϕ CD u^2 - (ϕ CD u^2)/(
   u + 1/α)) (Pi/CfD 1/(
    u (u/η + (2 Sqrt[u])/CfD)^(1/2)) + S/u)) /. datrul1;

qPlot1 = ListLogLogPlot[Table[
{t, NLInvDurb[qD1[u], u, t, NStep, imax, 0, 0.00000001]},
{t, ini, imax, istep}], PlotStyle -> Red, PlotLegends -> {α} /. datrul1, Joined -> True] // N

Running the code as is, will take a very long time. What I would like to do is have a varying time step on a logrithmic scale. So say 10 points for each log cycle. i.e. 0.001,0.002,...,0.009,0.01,0.02,...,0.09,0.1,...

Any ideas on how to handle this are much appreciated.

Just as a side, I have been using a simpler function to test the program, so you don't have to spend so much time with the long equation:

F[u_] = 1/(u^2 + 1);
f[t_] = Sin[t];

DurbTab = ListPlot[Table[
{t, NLInvDurb[F[u], u, t, 1000, 3 Pi, 0, 0.00000001]},
{t, 0, 3 Pi, 1/10}], PlotStyle -> Black, PlotRange -> Automatic, PlotLegends -> {"Approximation"} // N];
RealTab = Plot[{f[r]}, {r, 0, 3 Pi}, PlotStyle -> Blue, PlotLegends -> {"Real Solution"}];
Show[{DurbTab, RealTab}, GridLines -> Automatic]

Thanks

$\endgroup$
  • $\begingroup$ Welcome to Mathematica.SE! I suggest the following: 0) Browse the common pitfalls question 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Read the faq! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! $\endgroup$ – Dr. belisarius Mar 3 '16 at 19:37
  • $\begingroup$ Perhaps, use integer index i, and define t in the first argument of Table as 10^(i/10-3). $\endgroup$ – bbgodfrey Mar 3 '16 at 20:20
2
$\begingroup$

You can generate a power law list of indices with

ini = 0.001;
imax = 1000;
indices = 
  Append[Flatten @ Transpose @ Table[k PowerRange[ini, imax - 1], {k, 9}], 
  imax]
{0.001, 0.002, 0.003, 0.004, 0.005, 0.006, 0.007, 0.008, 0.009, 
 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 
 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 
 1., 2., 3., 4., 5., 6., 7., 8., 9., 
 10., 20., 30., 40., 50., 60., 70., 80., 90., 
 100., 200., 300., 400., 500., 600., 700., 800., 900., 
 1000}

and then index your table with

{t, indices}
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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