2
$\begingroup$

I need to plot the below equations:

$W^{+}(x,\tau,\alpha)]=-\frac{2m}{\pi\alpha}\exp[{\alpha x}]\times\int_{0}^{\infty}\frac{G(\omega,\tau_c)}{1+\omega^2}\cos(2\tau_c\omega\sqrt{1+\omega^2}-2\tan^{-1}(\omega) + \alpha x\omega)\text{d}\omega + 2mx$

and

$W^{-}(x,\tau,\alpha)]=-\frac{2m}{\pi\alpha}\exp[{\alpha x}]\times\int_{0}^{\infty}\frac{G(\omega,\tau_c)}{1+\omega^2}\cos(2\tau_c\omega\sqrt{1+\omega^2}-2\tan^{-1}(\omega) + \alpha x\omega)\text{d}\omega $

Here, $\tau_c=\tau\alpha^3$; $G(\omega,\tau_c)=\exp[\tau_c\sqrt{1+\omega^2}(1-\omega^2)]$ and $m=0.2$.

For $\alpha=0$ the profile should be symmetric, and for increasing $\alpha$ the profile becomes asymmetric across the x-axis. Here is my naive attempt

 G[\[Omega]_,\[Tau]C_] := Exp[\[Tau]C Sqrt[1+\[Omega]^2](1-\[Omega]^2)]   
  position[x_] := -2(mSlope/\[Alpha] \[Pi])Exp[\[Alpha] x] NIntegrate[G[\[Omega],\[Tau]C] Cos[2 \[Tau]C \[Omega] Sqrt[1+\[Omega]^2]-2ArcTan[\[Omega]]+\[Alpha] x \[Omega] ]/(1+\[Omega]^2),{\[Omega],0,\[Infinity]}] + 2 mSlope x
    positionLeft[x_] := -2(mSlope/\[Alpha] \[Pi])Exp[\[Alpha] x] NIntegrate[G[\[Omega],\[Tau]C] Cos[2 \[Tau]C \[Omega] Sqrt[1+\[Omega]^2]-2ArcTan[\[Omega]]+\[Alpha] x \[Omega] ]/(1+\[Omega]^2),{\[Omega],0,\[Infinity]}] 
\[Tau]C = 2; \[Alpha] = 0.001;\[mSlope] = 0.2
Plot[positionLeft[k], {k, 0.0, -10}]
Plot[position[k], {k, 0.0, 10}]

Output

enter image description here enter image description here

At $x=0$, the above plot does not show discontinuity in the profile slopes.

Can anyone help out on how to plot these type of equations and what type of NInegrate scheme one should apply?

P.S: The above figure and equations are from Eq. 5 from this article.

Thanks

$\endgroup$
4
  • $\begingroup$ Could you give a reference for where you got the plot in your question? $\endgroup$ May 19, 2020 at 11:51
  • $\begingroup$ @J.M. I have included the Ref. $\endgroup$
    – newstudent
    May 19, 2020 at 11:57
  • $\begingroup$ @newstudent What did you try so far yourself in Mathematica? $\endgroup$ May 19, 2020 at 14:42
  • $\begingroup$ @RolfMertig I have edited the post with my attempt and obtained output $\endgroup$
    – newstudent
    May 19, 2020 at 15:21

1 Answer 1

4
$\begingroup$

In your equations (mSlope/α π) should be (mSlope/(α π)). Your value of τC = 2 appears questionable since τC == τ*α^3 and with α = 10^-3, τ would be very large. I have used τC = 1*^-6

Clear["Global`*"]

G[ω_, τC_] := Exp[τC Sqrt[1 + ω^2] (1 - ω^2)]

position[x_] := Module[{xx = SetPrecision[x, 20]},
  -2 (mSlope/(α π)) Exp[α xx] NIntegrate[
     G[ω, τC] Cos[
        2 τC ω Sqrt[1 + ω^2] - 
         2 ArcTan[ω] + α xx ω]/(1 + ω^2), {\
ω, 0, ∞}, WorkingPrecision -> 20] + 2 mSlope xx]

positionLeft[x_] := Module[{xx = SetPrecision[x, 20]},
  -2 (mSlope/(α π)) Exp[α xx] NIntegrate[
    G[ω, τC] Cos[
       2 τC ω Sqrt[1 + ω^2] - 
        2 ArcTan[ω] + α xx ω]/(1 + ω^2), \
{ω, 0, ∞}, WorkingPrecision -> 20]]

W[x_] := Piecewise[{{position[x], x >= 0}, {positionLeft[x], x < 0}}]

τC = 1*^-6; α = 10^-3; mSlope = 1/5;

Plot[W[k]/mSlope, {k, -60, 60}, PlotRange -> All]

enter image description here

$\endgroup$
1
  • $\begingroup$ Thank you very much ! $\endgroup$
    – newstudent
    May 20, 2020 at 7:04

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.