Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Laplace Transform and complex contours

  1. Nov 2, 2005 #1

    saltydog

    User Avatar
    Science Advisor
    Homework Helper

    In the HW section, the following equation was proposed:

    [tex]ty(t)=\int_0^t \tau^{1-\alpha}y(t-\tau)d\tau;\quad \int_0^\infty y(x)dx=1[/tex]

    Link to HW thread: https://www.physicsforums.com/showthread.php?t=97562

    Using Laplace Transforms:

    [tex]Y(p)=\mathcal{L}\left\{y(t)\right\}[/tex]

    and noting the RHS is the inverse transform of a convolution, we obtain the IVP:

    [tex]-Y^{'}=Y\frac{\Gamma(\alpha)}{p^\alpha};\quad Y(0)=\int_0^\infty e^{-0t}y(t)dt=1[/tex]

    solving for [itex]\alpha \ne 1[/itex]:

    [tex]Y(p)=Exp\left[-\frac{\Gamma(\alpha)}{1-\alpha}s^{1-\alpha}\right][/tex]

    I (Mathematica) can invert this only for [itex]\alpha=1/2[/itex]:

    [tex]y(t)=\frac{e^{-\pi/t}}{t^{3/2}}[/tex]

    I'd like to know how to invert it for other values of alpha. My only recourse is to evaluate the complex contour for the inverse transform:

    [tex]\frac{1}{2\pi i}\int_{\gamma-i\infty}^{\gamma+i\infty} Y(p)e^{pt}dp[/tex]

    I'm not too good at this. I know I need to learn more about Complex Analysis and contours specifically, but for the time being, can someone help me through the integral for alpha=1/2:

    [tex]\frac{1}{2\pi i}\int_{\gamma-i\infty}^{\gamma+i\infty} \frac{e^{st}}{e^{2\sqrt{\pi s}}}ds=\frac{e^{-\pi/t}}{t^{3/2}}[/tex]
     
    Last edited: Nov 2, 2005
  2. jcsd
  3. Nov 2, 2005 #2

    Tide

    User Avatar
    Science Advisor
    Homework Helper

    Briefly: Make a change of variables using [itex]s = \sigma^2[/itex], complete the square in the exponential, scale out the t (which gives you the scaling of [itex]t^{-3/2}[/itex] immediately) and note that the change of variables requires your new contour to asymptote to the lines [itex]e^{i \pi / 4}[/itex] and [itex]e^{-i \pi / 4}[/itex]. It should be apparent that you're dealing with the Fresnel integral functions with the usual results.
     
  4. Nov 2, 2005 #3

    saltydog

    User Avatar
    Science Advisor
    Homework Helper

    Thanks Tide. I'll work with it. May take a while. :smile:
     
  5. Nov 2, 2005 #4

    saltydog

    User Avatar
    Science Advisor
    Homework Helper

    Tide, it's not happening for me. Doing as you suggest, (completing the square), I get:

    [tex]2e^{-\pi/t^2}\int e^{t(u-\frac{\sqrt{\pi}}{t})^2}udu[/tex]

    I don't see where to go from there and I really don't understand how the limits are changed as you suggest.

    But I know where my problem lies and I know how to correct it: I don't have a sufficient understanding of Complex Analysis to even ask the question, and in order to approach this problem I would first have to work on other more simple problems. That's Ok. I have the time.

    I did however go to the Library today and searched several tables of inverse transforms. It turns out that only for the case of alpha=1/2 is an inverse transform published. I suspect then even when I do master the contour integral, a closed-form solution for other values of alpha will NOT be attainable and thus I will have to resort to numerical methods of computing the inverse transforms. I know absolutely nothing about such numerical methods but I'm interested in learning about them and arriving at a complete solution of this problem in terms of alpha.

    I also purchased a book on Complex Analysis and will take my own advice in here: Start on page 1.:smile:
     
  6. Nov 3, 2005 #5

    Tide

    User Avatar
    Science Advisor
    Homework Helper

    salty,

    I just went through the transformations and I get

    [tex]\frac {2 e^{-\pi / t}}{t^{3/2}} \int e^{\left( \sigma - \sqrt {\pi / t}\right)^2} \sigma d\sigma[/tex]

    The original contour for the inverse Laplace transform effectively goes from [itex]R e^{-i\pi/2}[/itex] to [itex]R e^{i\pi/2}[/itex] with [itex]R \rightarrow \infty[/itex]. Since [itex]\sigma = \sqrt s[/itex] the new limits will be from [itex]\sqrt R e^{-i\pi/4}[/itex] to [itex]\sqrt R e^{i\pi/4}[/itex]. Looking at it more carefully shows one more change of variables is needed (a translation) to simplify the argument of the exponential.

    When you're all one with that, you should have the integral I suggested in my previous post. For now you just have to go back and recheck your transformations.
     
    Last edited: Nov 3, 2005
  7. Nov 3, 2005 #6

    Tide

    User Avatar
    Science Advisor
    Homework Helper

    P.S. I didn't include the factor [itex]\frac {1}{2\pi i}[/itex] in my previous post.
     
  8. Nov 3, 2005 #7

    saltydog

    User Avatar
    Science Advisor
    Homework Helper

    Tide, I cannot arrive at the answer you propose and when I check the integral with real values, your substitution does not lead to the numerical answer:

    Let's just for now consider a real integral and let t=3:

    [tex]\int_0^1 \frac{e^{st}}{e^{2\sqrt{\pi s}}}ds[/tex]

    when I numerically integrate this I obtain:

    [tex]\int_0^1 \frac{e^{st}}{e^{2\sqrt{\pi s}}}ds\approx 0.428599[/tex]

    When I let:

    [tex]\sigma^2=s[/tex]

    then the limits in [itex]\sigma[/itex] are the same.

    Completing the square as you suggest, I arrive at:

    [tex]2e^{-\pi/t}\int_0^1 e^{t(\sigma-\sqrt{\pi}/t)^2}\sigma d\sigma[/tex]

    numerically integrating that I obtain the same result as above. However, when I numerically integrate your solution:

    [tex]\frac {2 e^{-\pi / t}}{t^{3/2}} \int_0^1 e^{\left( \sigma - \sqrt {\pi / t}\right)^2} \sigma d\sigma[/tex]

    I obtain 0.0830375. I suspect I'm making a dumb mistake and can handle you correcting me.
     
    Last edited: Nov 3, 2005
  9. Nov 3, 2005 #8

    Tide

    User Avatar
    Science Advisor
    Homework Helper

    salty,

    There are two problems. First, I didn't do the transformations quite right so the [itex]t^(3/2)[/itex] should just be [itex]t[/itex]:

    [tex]\frac {2 e^{-\pi / t}}{t} \int e^{\left( \sigma - \sqrt {\pi / t}\right)^2} \sigma d\sigma[/tex]

    Also, if you want do a numerical comparison you have to set the limits properly. Remember, there was a scaling that eliminated the factor of t from the exponential which looks like [itex]\sigma = \sqrt {t} \sigma'[/itex] so if you set t = 3 then the upper limit in the transformed integral will be [itex]\sqrt {3}[/itex].
     
  10. Nov 3, 2005 #9

    saltydog

    User Avatar
    Science Advisor
    Homework Helper

    Tide,

    As you requested, I've achieved oneness with the integral. Here it is with your substitutions in case others are following this:

    Above, after letting:

    [tex]\sigma^2=s[/tex]

    and completing the square, we obtain:

    [tex]2e^{-\pi/t}\int e^{t(\sigma-\sqrt{\pi}/t)^2}\sigma d\sigma[/tex]

    Now, in order to remove the t in the exponent, we let:

    [tex]v=\sqrt{t}\sigma[/tex]

    so that:

    [tex]dv=\sqrt{t}d\sigma,\quad \sigma=\frac{v}{\sqrt{t}},\quad d\sigma=\frac{dv}{\sqrt{t}}[/tex]

    Substituting this scalling factor into the exponent:

    [tex]
    \begin{align*}
    t\left[\frac{v^2}{t}-\frac{2v\sqrt{\pi}}{t\sqrt{t}}+\frac{\pi}{t^2}\right]&=
    v^2-2v\sqrt{\pi/t}+\frac{\pi}{t} \\
    &=(v-\sqrt{\pi/t})^2
    \end{align}
    [/tex]

    substituting this into the integral:

    [tex]2e^{-\pi/t}\int e^{(v-\sqrt{\pi/t})^2}\left(\frac{v}{\sqrt{t}}\right)\frac{dv}{\sqrt{t}}[/tex]

    Simplifying:

    [tex]\frac{2e^{-\pi/t}}{t}\int e^{(v-\sqrt{\pi/t})^2}vdv;\quad v=\sqrt{t}\sigma;\quad \sigma^2=s[/tex]

    I'll start working with the complex analysis part . . .
     
    Last edited: Nov 3, 2005
  11. Nov 4, 2005 #10

    saltydog

    User Avatar
    Science Advisor
    Homework Helper

    Tide, anyone else too you know, can you please review the analysis below and let me know if I understand it so far. Thanks!

    The original integral was:

    [tex]\int_{\gamma-i\infty}^{\gamma+i\infty} \frac{e^{st}}{e^{2\sqrt{\pi s}}}ds[/tex]

    This can be re-written as:

    [tex]\mathop\lim\limits_{y\to\infty}\int_{\gamma-yi}^{\gamma+yi} \frac{e^{st}}{e^{2\sqrt{\pi s}}}ds[/tex]

    Gamma is taken to be ANY positive constant such that the contour is to the right of any singularities of the integrand. I believe the integrand is entire so we can just take gamma to be 1. Converting to polar coordinates, the upper and lower limits become:

    [tex]s_1=1-yi=Re^{-i\theta}[/tex]

    [tex]s_2=1+yi=Re^{i\theta}[/tex]

    Now, as y goes to infinity, theta approaches pi/2 and R goes to infinity. In the first transformation we had:

    [tex]\sigma^2=s[/tex]

    Thus in terms of sigma, the limits become:

    [tex]\sigma_1=\sqrt{Re^{-\pi i/2}}=\sqrt{R}e^{-\pi i/4}[/tex]

    [tex]\sigma_2=\sqrt{Re^{\pi i/2}}=\sqrt{R}e^{\pi i/4}[/tex]

    with R approaching infinity.

    In the second transformation, we have:

    [tex]v=\sqrt{t}\sigma[/tex]

    In terms of v, then the limits become:

    [tex]v_1=\sqrt{Rt}e^{-\pi i/4}[/tex]

    [tex]v_2=\sqrt{Rt}e^{\pi i/4}[/tex]

    Thus arriving at the final form of the inverse transform:

    [tex]y(t)=\frac{1}{2\pi i}\frac{2e^{-\pi/t}}{t}\mathop\lim\limits_{R\to\infty}\int_{\sqrt{Rt}e^{-\pi i/4}}^{\sqrt{Rt}e^{\pi i/4}} e^{(v-\sqrt{\pi/t})}^2vdv[/tex]

    I've already done some preliminary calculations of the integral and the results do indeed tend to the desired result:

    [tex]y(t)=\frac{e^{-\pi/t}}{t^{3/2}}[/tex]

    However I cannot yet provide a proof and plan to work on this next.:smile:
     
  12. Nov 4, 2005 #11

    saltydog

    User Avatar
    Science Advisor
    Homework Helper

    For the following:

    [tex]y(t)=\frac{1}{2\pi i}\frac{2e^{-\pi/t}}{t}\mathop\lim\limits_{R\to\infty}\int_{\sqrt{RT}e^{-\pi i/4}}^{\sqrt{RT}e^{\pi i/4}} e^{(v-\sqrt{\pi/t})}^2vdv[/tex]

    First consider:

    [tex]\int e^{(x-k)^2} xdx[/tex]

    Letting:

    [tex]u=x,\quad dv=e^{(x-k)^2}dx[/tex]

    so that:

    [tex]du=dx,\quad v=\frac{\sqrt{\pi}}{2}\text{Erfi}(x-k)[/tex]

    where:

    [tex]\text{Erfi}(z)=\frac{\text{Erf}(iz)}{i}[/tex]

    and therefore:

    [tex]\int e^{(x-k)}^2 xdx=x\frac{\sqrt{\pi}}{2}\text{Erfi}(x-k)-\frac{\sqrt{\pi}}{2}\int \text{Erfi}(x-k)dx[/tex]

    Now:

    [tex]\int \text{Erfi}(x-k)dx=\frac{1}{\sqrt{\pi}}e^{(x-k)^2}-k\text{Erfi}(x-k)+x\text{Erfi}(x-k)[/tex]

    and thus:

    [tex]\int e^{(x-k)^2}xdx=\frac{\sqrt{\pi}}{2}\left[\frac{e^{(x-k)^2}}{\sqrt{\pi}}+k\text{Erfi}(x-k)\right][/tex]

    Therefore if we let:

    [tex]k=\sqrt{\pi/t}[/tex]

    we obtain:

    [tex]\int e^{(v-\sqrt{\pi/t})^2}vdv=\frac{\sqrt{\pi}}{2}\left[\frac{e^{(v-\sqrt{\pi/t})^2}}{\sqrt{\pi}}+\sqrt{\pi/t}\text{Erfi}(v-\sqrt{\pi/t})\right][/tex]

    This then reduces the inverse transform to the following expression:

    [tex]y(t)=\frac{e^{-\pi/t}}{2\pi i t}\left[\mathop\lim\limits_{R\to\infty}\left(e^{(v-\sqrt{\pi/t})^2}+\frac{\pi}{t^{1/2}}\text{Erfi}(v-\sqrt{\pi/t})\right)_{\sqrt{Rt}e^{-\pi i/4}}^{\sqrt{Rt}e^{\pi i/4}}\right][/tex]

    Edit: Oh yea: Equal rights for special functions.
     
    Last edited: Nov 4, 2005
  13. Nov 5, 2005 #12

    saltydog

    User Avatar
    Science Advisor
    Homework Helper

    I wish to analyze:

    [tex]y(t)=\frac{e^{-\pi/t}}{2\pi i t}\left[\mathop\lim\limits_{R\to\infty}\left(e^{(v-\sqrt{\pi/t})^2}+\frac{\pi}{t^{1/2}}\text{Erfi}(v-\sqrt{\pi/t})\right)_{\sqrt{Rt}e^{-\pi i/4}}^{\sqrt{Rt}e^{\pi i/4}}\right][/tex]

    in two parts.

    First consider the limit of the first term:

    [tex]\mathop\lim\limits_{R\to\infty}\left[e^{(v-\sqrt{\pi/t})^2}\right]_{\sqrt{RT}e^{-\pi i/4}}^{\sqrt{Rt}e^{\pi i/4}}[/tex]

    Converting the upper limit to polar form:

    [tex]
    \begin{align*}
    \sqrt{Rt}e^{\pi i/4}&=\sqrt{Rt}\left(Cos(\pi/4)+iSin(\pi/4)\right) \\
    &=\sqrt{\frac{Rt}{2}}+i\sqrt{\frac{Rt}{2}}
    \end{align}
    [/tex]

    so that:

    [tex]
    \begin{align*}
    (v-\sqrt{\pi/t})^2&=\left[\left(\sqrt{Rt/2}-\sqrt{\pi/2}\right)+i\sqrt{Rt/2}\right]^2 \\
    &=\left(-2\sqrt{Rt/2}+\pi/t\right)+i\left(Rt-2\sqrt{R\pi/2}\right)
    \end{align}
    [/tex]

    The upper limit then becomes:

    [tex]\mathop\lim\limits_{R\to\infty} Exp\left[\left(-2\sqrt{R\pi/2}+\pi/t\right)+i\left(Rt-2\sqrt{R\pi/2}\right)\right][/tex]

    Note the real part of the expression is negative and as R goes to infinity, the coefficient on the Euler expansion goes to zero and thus the limit is zero.

    The same analysis can be made for the lower limit. This limit is therefore zero. We can now write the inverse transform as:

    [tex]y(t)=\frac{e^{-\pi/t}}{2\pi i t}\mathop\lim\limits_{R\to\infty}\left\{\frac{\pi}{t^{1/2}}\text{Erfi}\left(v-\sqrt{\pi/t}\right)_{\sqrt{Rt}e^{-\pi i/4}}^{\sqrt{Rt}e^{\pi i/4}}\right\}[/tex]

    or:

    [tex]y(t)=\frac{e^{-\pi/t}}{2 i t^{3/2}}\mathop\lim\limits_{R\to\infty}\left\{\text{Erfi}\left(v-\sqrt{\pi/t}\right)_{\sqrt{Rt}e^{-\pi i/4}}^{\sqrt{Rt}e^{\pi i/4}}\right\}[/tex]
     
    Last edited: Nov 5, 2005
  14. Nov 5, 2005 #13

    saltydog

    User Avatar
    Science Advisor
    Homework Helper

    I wish to finally evaluate:

    [tex]
    y(t)=\frac{e^{-\pi/t}}{2 i t^{3/2}}\mathop\lim\limits_{R\to\infty}\left\{\text{Erfi}\left(v-\sqrt{\pi/t}\right)_{\sqrt{Rt}e^{-\pi i/4}}^{\sqrt{Rt}e^{\pi i/4}}\right\}\tag{1}
    [/tex]

    First note:

    [tex]
    \text{Erfi}(z)=\frac{\text{Erf(iz)}}{i}=\frac{1}{i}\frac{2}{\sqrt{\pi}}\int_0^{iz} e^{-w^2}dw
    [/tex]

    Using this relation, I can re-write the above limit as:

    [tex]
    \frac{2}{i\sqrt{\pi}}\left\{\mathop\lim\limits_{R\to\infty}\int_0^{i(\sqrt{Rt}e^{\pi i/4})} e^{-w^2}dw-\mathop\lim\limits_{R\to\infty}\int_0^{i(\sqrt{Rt}e^{-\pi i/4})} e^{-w^2}dw\right\}
    [/tex]

    Using Euler's relation, I convert the upper limits to real and imaginary parts and multiply by i:

    [tex]
    \begin{align*}
    i\sqrt{Rt}e^{\pi i/4}&=i\left[\left(\sqrt{Rt/2}-\sqrt{\pi/2}\right)+i\sqrt{Rt/2}\right] \\
    &=-\sqrt{Rt/2}+i\left(\sqrt{Rt/2}-\sqrt{\pi/t}\right)
    \end{align}
    [/tex]

    Now, as R goes to infinity, this quantity goes to [itex]-\infty+i\infty[/itex].

    And the lower limit:

    [tex]
    \begin{align*}
    i\sqrt{Rt}e^{-\pi i/4}&=i\left[\sqrt{Rt/2}+\left(\sqrt{Rt/2}-\sqrt{\pi/2}\right)\right] \\
    &=\sqrt{Rt/2}+i\left(\sqrt{Rt/2}-\sqrt{\pi/2}\right)
    \end{align}
    [/tex]

    as R goes to infinity, this limit goes to [itex]\infty+i\infty[/itex].

    Thus the integrals can be rewritten as:

    [tex]
    \frac{1}{i}\frac{2}{\sqrt{\pi}}\int_0^{-\infty+i\infty} e^{-w^2}dw-
    \frac{1}{i}\frac{2}{\sqrt{\pi}}\int_0^{\infty+i\infty} e^{-w^2}dw
    [/tex]

    I choose to evaluate these line integrals along a parametric path using the following relation from Complex Analysis:

    [tex]
    \int_C f(z)dz=\int_c udx-vdy+i\int_c vdx+udy
    [/tex]

    Thus the first step is to express the integrand in terms of a real and imaginary part. Relying on Euler's formula:

    [tex]
    \begin{align*}
    e^{-w^2}&=e^{-(x+yi)^2} \\
    &=e^{-(x^2-y^2+2xyi)} \\
    &=e^{-(x^2-y^2)}Cos(2xy)+ie^{-(x^2-y^2)}(-Sin(2xy))\\
    &=u+iv
    \end{align}
    [/tex]

    For the first integral, the path goes from the origin to the point [itex]-\infty+i\infty[/itex]

    Thus if I let:

    [tex]x=-t,\quad dx=-dt[/tex]

    [tex]y=t,\quad dy=dt[/tex]

    Making this substitution, note that the exponent term is 1 and we have for the first integral (with t going from 0 to infinity):

    [tex]
    \begin{align*}
    \int_0^{-\infty+i\infty} e^{-w^2}dw&=\left(-\int_0^\infty Cos(-2t^2)dt+\int_0^\infty Sin(-2t^2)dt\right) \\
    &+i\left(\int_0^\infty Cos(-2t^2)dt+\int_0^\infty Sin(-2t^2)dt\right) \\
    &=\left(-\frac{\sqrt{\pi}}{4}-\frac{\sqrt{\pi}}{4}\right)+i\left(\frac{\sqrt{\pi}}{4}-\frac{\sqrt{\pi}}{4}\right) \\
    &=-\frac{\sqrt{\pi}}{2}
    \end{align}
    [/tex]


    The same analysis can be made (with the path going from 0 to [itex]\infty+i\infty[/itex] for the second integral with the results that the line integral in that case is:

    [tex]
    \frac{\sqrt{\pi}}{2}[/tex]

    Substituting these values into (1), I obtain:

    [tex]
    \frac{1}{i}\frac{2}{\sqrt{\pi}}\int_0^{-\infty+i\infty} e^{-w^2}dw=\frac{1}{i}\frac{2}{\sqrt{\pi}}(-\frac{\sqrt{\pi}}{2})=-\frac{1}{i}\frac{i}{i}=i
    [/tex]

    and:

    [tex]
    \frac{1}{i}\frac{2}{\sqrt{\pi}}\int_0^{\infty+i\infty} e^{-w^2}dw=\frac{1}{i}\frac{2}{\sqrt{\pi}}(\frac{\sqrt{\pi}}{2})=\frac{1}{i}\frac{i}{i}=-i
    [/tex]

    So that the limit is [itex]2i[/itex].


    Substituting this into the last expression of the inverse transform, we obtain:

    [tex]
    y(x)=\frac{e^{-\pi/t}}{2i t^{3/2}}(2i)=\frac{e^{-\pi/t}}{t^{3/2}}[/tex]

    Thanks for helping me Tide. Couldn't have done this without your help.:smile:
     
  15. Nov 5, 2005 #14

    Tide

    User Avatar
    Science Advisor
    Homework Helper

    Way to go, Salty!
     
  16. Nov 6, 2005 #15

    saltydog

    User Avatar
    Science Advisor
    Homework Helper

    So I wish to return to the subject of differential (integral) equations and the numerical evaluation of inverse Laplace Transforms. While going though the analysis above, I realized I can numerically evaluate an integral of a complex variable in the same (similar) way as for a real variable! Take for example the contour expression for the inverse transform:

    [tex]y(x)=\frac{1}{2\pi i}\int_{1-i\infty}^{1+i\infty}\frac{e^{st}}{e^{2\sqrt{\pi s}}} ds[/tex]

    For the moment, let's consider the case of t=1. In that case, this integral is EXACTLY equal to:

    [tex]y(1)=e^{-\pi}\approx 0.0432139[/tex]

    Now, for starters, suppose I just calculate a ROUGH Reimann sum for this integral:

    [tex]\frac{1}{2\pi i}\sum_1^{1000} \frac{e^{st}}{e^{2\sqrt{\pi s}}}\Delta c[/tex]

    along the straight-line contour from the points [itex]1-10i[/itex] to [itex]1+10i[/itex].

    This is my Mathematica code:

    Code (Text):
    sum = 0;
    tval = 1;
    imax = 10;
    cval = 1 - i max;
    nmax = 1000;
    deltai = 2 imax/nmax i;
    For[n = 1, n <= nmax, n++,
      cval += deltai;
    [tex]sum+=N\left[\left(\frac{e^{\text{cval tval}}}{e^{2\sqrt{\pi\text{cval}}}}\right) \text{deltai}\left][/tex]
    ];
    [tex]sum=\frac{1}{2\pi i} sum[/tex]
    The results are:

    [tex]0.0434844+1.31273X10^{-6} i[/tex]

    Granted there is a slight complex component but the real part agrees with the actual value to 3 decimal places and this is only a rough estimate.

    My suspicion is that this will be the case ONLY for alpha=1/2 and that any other value of alpha in the original integral equation will result in a significant complex part meaning that the only REAL solution to the integral equation is the one with alpha=1/2.

    Anyone have any thoughts about this?

    Tide, you still with me?
     
    Last edited: Nov 6, 2005
  17. Nov 6, 2005 #16

    Tide

    User Avatar
    Science Advisor
    Homework Helper

    If real solutions exist to the original integral equation and if the Laplace transforms exist then the inversion must lead to a real result. A numerical approximation may lead to an imaginary part but there are two issues with a numerical approximation: (a) finite step size and (b) infinite domain. It may take that entire "trip" to infinity to eliminate the small imaginary part but it will not be there in the end! :)
     
  18. Nov 6, 2005 #17

    saltydog

    User Avatar
    Science Advisor
    Homework Helper

    Numerical calculations

    Hey Tide, I find the following results very encouraging: The first plot is the solution of the integral equation:

    [tex]y(t)=\frac{e^{-\pi/t}}{t^{3/2}}[/tex]

    The second plot is a numerical evaluation of the inverse transform using a simple Riemann sum over the path [itex]1-20i[/itex] to [itex]1+20i[/tex] sampling 1000 points.

    The third plot is a superposition of the two. Not bad for a rough approximation.:smile:

    Edit: I neglected the "small" residual imaginary component of the numerical calculations.
     

    Attached Files:

    Last edited: Nov 6, 2005
  19. Nov 6, 2005 #18

    Tide

    User Avatar
    Science Advisor
    Homework Helper

    That's pretty good!

    If you want to take it a little further, you could do some asymptotic analysis on the integral. When t is large, most of the contribution to the integral (say your first eq. in #11) is from the region about a "saddle point" (e.g. [itex]\int e^{\phi(s, t)} ds[/itex]) has such a point where [itex]d \phi(s, t) /ds = 0[/itex] and you integrate along the path of steepest descent through the saddle point.

    In principle, the asymptotic expansion of the integral will overlap your numerical evaluation over some range of large t values. I haven't looked but you should be able to find a lot of references to the "method of steepest descent" for details.
     
    Last edited: Nov 6, 2005
  20. Nov 7, 2005 #19

    saltydog

    User Avatar
    Science Advisor
    Homework Helper

    Thanks Tide. I'll check that out. Here are the plots for alpha =0.5 to alpha=0.1 with the largest one being alpha=0.1, then 0.2, 0.3, 0.4, and the smallest, 0.5.

    I've yet to show that these other plots satisfy the integral equation however and plan to work on a back-substitution method to check them (within some level of accuracy).

    I'm sure other numerical methods exists for evaluating integral equations, however, to pursue such would have denied me the priviledge (and learning experience) of working on this very interesting approach with Laplace Transforms.:smile:
     

    Attached Files:

  21. Nov 8, 2005 #20

    saltydog

    User Avatar
    Science Advisor
    Homework Helper

    I perfected a method of back-substitution of the numeric data into the integral equation in an effort to numerically check the solution for alpha=0.1. The inverse transform is:

    [tex]y(t)=\frac{1}{2\pi i}\mathop\lim\limits_{a\to\infty}\int_{1-ai}^{1+ai} e^{zt}\text{Exp}\left[-\frac{\Gamma[0.1]}{1-0.1}z^{1-0.1}\right]dz[/tex]

    This can be numerically integrated directly until the results reach a desired level of accuracy. The integral then is evaluated for a range of t. The results of these calculations is shown in plot 1. Mathematica has a nice feature which can represents the data points as a "Interpolation function":

    [tex]\text{alpha01 = Interpolation[pointlist]}[/tex]

    alpha01 is now treated as a function.

    The next step is to numerically evaluate the RHS of the original integral equation:

    [tex]\int_0^t \tau^{\alpha-1}\text{alpha01}(t-\tau)d\tau[/tex]

    or:

    [tex]rhs[t]=\int_0^t (t-z)^{\alpha-1}\text{alpha01}(z)dz[/tex]

    The LHS of the integral equation is:

    [tex]lhs[t]=\text{alpha01[t]t}[/tex]

    Subtracting:

    [tex]\text{dif[t]=lhs[t]-rhs[t]}[/tex]

    The results of this difference is shown in the second plot showing that the difference is no larger than 10^-6.

    This then gives me some confidence that the numerical results for alpha=0.1 do indeed satisfy the integral equation.
     

    Attached Files:

    Last edited: Nov 8, 2005
  22. Nov 10, 2005 #21

    saltydog

    User Avatar
    Science Advisor
    Homework Helper

    For the record I wish to correct a statement I made above:

    The integrand has a branch point (singularity) at the origin (caused by the square root of z). Pretty sure anyway. Anyone reading this thread, please forgive my unfamiliarity with the subject. This is all very new for me.
    With this in mind, I'd like to evaluate the inverse transform using Cauchy's theorem:

    [tex]\oint_c \frac{e^{zt}}{e^{2\sqrt{\pi z}}}dz=0[/tex]

    along the contour shown in the first plot below. The second plot is an expanded view of the center contour from D to E. You guys mind if I continue working on this here even though the contour integration is not "directly" relevant to differential equations? The integrals are tough to evaluate and I'll likely need help with them. I've already done some preliminary numerical calculations with all 5 integrals and the results do indeed tend to the desired result.
     

    Attached Files:

    Last edited: Nov 10, 2005
  23. Nov 11, 2005 #22

    saltydog

    User Avatar
    Science Advisor
    Homework Helper

    Before I begin evaluating this contour integral I'd like to justify the path with my limited understanding of the subject:

    I need to avoid going through the singularity of a branch cut created by the z^{1/2} term (don't quite understand this part). This cut can be anywhere in the complex plane. Since the definition of the inverse transform requires the path of integration to be right of any singularities, for simplicity I define the cut to be the negative x-axis. Thus to make a closed contour, we go around both the branch cut and the branch point at the origin, hense the keyhole contour. For those unfamiliar with the subject, the actual integration paths are infinitessimally close to both the negative x-axis and the origin. We simply expand the path in the contour diagram so that we can more easily visualize it.

    Since the branch cut is at the polar coordinate [itex]\theta=\pi[/itex], in any integrations involving polar coordinates, we must not go through the value [itex]\theta=\pi[/itex]. Thus in the calculations of the integrals, theta will take on values between 0 and pi and between -pi and pi.

    In order to facillitate the evaluation of the inverse transform, we allow the diameter (R in the diagram) of the contour to approach infinity and evaluate the integral at that limit.

    Thus, using Cauchy's Integral Theorem and the contour defined above, I integrate along a counter-clockwise direction and arrive at the following expression:

    [tex]
    \begin{align*}
    \oint_c f(z)dz&=\int_{AB} f(z)dz+\int_{BC} f(z)dz+\int_{CD} f(z)dz\\
    &+ \int_{DE} f(z)dz+\int_{EF} f(z)dz+\int_{FA} f(z)dz=0
    \end{align}
    [/tex]

    Note however at the limit [itex]R\to\infty[/itex], that

    [tex]\int_{AB} f(z)dz=\int_{\gamma-i\infty}^{\gamma+i\infty} f(z)dz[/tex]

    This is in fact the integral component of the transform. Thus I can write:

    [tex]
    \begin{align*}
    \int_{\gamma-i\infty}^{\gamma+i\infty} f(z)dz&=-\left(\int_{BC} f(z)dz+\int_{CD} f(z)dz+\int_{DE} f(z)dz\right \\
    &+\left\int_{EF} f(z)dz+\int_{FA} f(z)dz\right)
    \end{align}
    [/tex]

    Multiplying both sides by [itex]\frac{1}{2\pi i}[/itex], the inverse transform can be expressed as:

    [tex]
    \begin{align*}
    \frac{1}{2\pi i}\int_{\gamma-i\infty}^{\gamma+i\infty} f(z)dz&=-\frac{1}{2\pi i}\left(\int_{BC} f(z)dz+\int_{CD} f(z)dz+\int_{DE}
    f(z)dz\right \\
    &+ \left\int_{EF} f(z)dz+\int_{FA} f(z)dz\right)
    \end{align}
    [/tex]
     
    Last edited: Nov 11, 2005
  24. Nov 12, 2005 #23

    saltydog

    User Avatar
    Science Advisor
    Homework Helper

    Letting [itex]z=Re^{i\theta}[/itex] for the first and fourth integral, I obtain:

    [tex]\int_{BC}=\mathop\lim\limits_{R\to\infty}\int_{\pi/2}^{\pi}\frac{e^{tRe^{i\theta}}}{Exp[2\sqrt{R\pi}e^{i\theta/2}]} iRe^{i\theta}d\theta[/tex]

    [tex]\int_{FA}=\mathop\lim\limits_{R\to\infty}\int_{-\pi}^{-\pi/2}\frac{e^{tRe^{i\theta}}}{Exp[2\sqrt{R\pi}e^{i\theta/2}]} iRe^{i\theta}d\theta[/tex]

    For the center contour, let [itex]\phi[/itex] go from [itex]\pi[/itex] to[itex]-\pi[/itex] and let the modulus be [itex]\epsilon[/itex]. Thus letting:

    [tex]z=\epsilon e^{\phi i}[/tex]

    [tex]\int_{DE}=\mathop\lim\limits_{\epsilon\to 0}\int_{\pi}^{-\pi} \frac{e^{t\epsilon e^{i\phi}}}{Exp[2\sqrt{\pi\epsilon}e^{\phi i/2}]}\epsilon i e^{\phi i}d\phi[/tex]

    For the two straight-line segements, I'll use a common (apparently) substituion in complex integration:

    For the path C to D, let:

    [tex]z=-u=ue^{\pi i}[/tex]

    For the path E to F, let:

    [tex]z=-u=ue^{-\pi i}[/tex]

    Making these substitutions I obtain:

    [tex]\int_{CD}=\mathop\lim\limits_{R\to\infty,\epsilon\to 0}\int_R^{\epsilon}\frac{e^{-tu}}{e^{2\sqrt{\pi}i\sqrt{u}}} (-du}[/tex]

    [tex]\int_{EF}=\mathop\lim\limits_{R\to\infty,\epsilon\to 0}\int_{\epsilon}^{R}\frac{e^{-tu}}{e^{-2\sqrt{\pi}i\sqrt{u}}} (-du}[/tex]

    Putting this all together, I obtain for the inverse transform:

    [tex]
    \begin{align*}
    \frac{1}{2\pi i}\int_{\gamma-i\infty}^{\gamma+\infty}\frac{e^{zt}}{e^{2\sqrt{\pi z}}} dz&=
    -\frac{1}{2\pi i}\left(\mathop\lim\limits_{R\to\infty}\int_{\pi/2}^{\pi}\frac{e^{tRe^{i\theta}}}{Exp[2\sqrt{R\pi}e^{i\theta/2}]} iRe^{i\theta}d\theta \right\\ \\
    &+\mathop\lim\limits_{R\to\infty}\int_{-\pi}^{-\pi/2}\frac{e^{tRe^{i\theta}}}{Exp[2\sqrt{R\pi}e^{i\theta/2}]} iRe^{i\theta}d\theta \\ \\

    &+\mathop\lim\limits_{\epsilon\to 0}\int_{\pi}^{-\pi} \frac{e^{t\epsilon e^{i\phi}}}{Exp[2\sqrt{\pi\epsilon}e^{\phi i/2}]}\epsilon i e^{\phi i}d\phi \\ \\

    &+\mathop\lim\limits_{R\to\infty,\epsilon\to 0}\int_R^{\epsilon}\frac{e^{-tu}}{e^{2\sqrt{\pi}i\sqrt{u}}} (-du) \\ \\

    &+\left\mathop\lim\limits_{R\to\infty,\epsilon\to 0}\int_{\epsilon}^{R}\frac{e^{-tu}}{e^{-2\sqrt{\pi}i\sqrt{u}}}(-du)\right)


    \end{align}
    [/tex]

    Alright, so I know the first, second, and third go to zero but can't show that for the first and second (not sure about third), and the sum of the fourth and fifth go to the desired result and can kinda' get close to showing that. May take a while . . . can I use the HW section if I need help?

    Edit: Oh yea, I realize I'm bitting off more than I can chew with this. It's only a goal I wish to strive for; I'm currently working on simpler ones: Schaum's is a start.
     
    Last edited: Nov 12, 2005
  25. Nov 20, 2005 #24

    saltydog

    User Avatar
    Science Advisor
    Homework Helper

    Since I'm unable to evaluate the limits of the first two in terms of [itex]\theta[/itex], I'll rely on the analyticity of the integrand in the region of integration and it's indepenence of path.

    So I'll leave them as:

    [tex]\int_{BC} \frac{e^{zt}}{e^{\sqrt{\pi z}}}dz[/tex]

    [tex]\int_{FA} \frac{e^{zt}}{e^{\sqrt{\pi z}}}dz[/tex]

    and solve for an antiderivative exactly as I did above to obtain:

    [tex]\int_{BC} \frac{e^{zt}}{e^{\sqrt{\pi

    z}}}dz=\frac{2e^{-\pi/t}}{t}\mathop\lim\limits_{R\to\infty}\int_{\sqrt{Rt}e^{\pi

    i/4}}^{\sqrt{Rt}e^{\pi i/2}} e^{(v-\pi/t)^2}vdt[/tex]

    and:

    [tex]\int_{FA} \frac{e^{zt}}{e^{\sqrt{\pi

    z}}}dz=\frac{2e^{-\pi/t}}{t}\mathop\lim\limits_{R\to\infty}\int_{\sqrt{Rt}e^{-\pi

    i/2}}^{\sqrt{Rt}e^{-\pi i/4}} e^{(v-\pi/t)^2}vdt[/tex]

    The two integrals in z thus reduce to the same expression above (with different limits):

    [tex]\int_{BC} f(z)dz=\frac{2e^{-\pi/t}}{t}\left[\mathop\lim\limits_{R\to\infty}\left(e^{(v-\sqrt{\pi/t})^2}+\frac{\pi}{t^{1/2}}\text{Erf

    i}[v-\sqrt{\pi/t}]\right)_{\sqrt{Rt}e^{\pi i/4}}^{\sqrt{Rt}e^{\pi

    i/2}}\right][/tex]

    [tex]\int_{FA} f(z)dz=\frac{2e^{-\pi/t}}{t}\left[\mathop\lim\limits_{R\to\infty}\left(e^{(v-\sqrt{\pi/t})^2}+\frac{\pi}{t^{1/2}}\text{Erf

    i}[v-\sqrt{\pi/t}]\right)_{\sqrt{Rt}e^{-\pi i/2}}^{\sqrt{Rt}e^{-\pi

    i/4}}\right][/tex]

    Using the same analysis as above, these limits can be shown to approach zero.

    Thus I can reduce the expression of the inverse transform to:

    [tex]
    \begin{align*}
    \frac{1}{2\pi i}\int_{\gamma-i\infty}^{\gamma+\infty}\frac{e^{zt}}{e^{2\sqrt{\pi z}}} dz&=
    -\frac{1}{2\pi i}\left(\mathop\lim\limits_{\epsilon\to 0}\int_{\pi}^{-\pi} \frac{e^{t\epsilon e^{i\phi}}}{Exp[2\sqrt{\pi\epsilon}e^{\phi i/2}]}\epsilon i e^{\phi i}d\phi \\ \\
    &+\mathop\lim\limits_{R\to\infty,\epsilon\to 0}\int_R^{\epsilon}\frac{e^{-tu}}{e^{2\sqrt{\pi}i\sqrt{u}}} (-du) \\ \\
    &+\left\mathop\lim\limits_{R\to\infty,\epsilon\to 0}\int_{\epsilon}^{R}\frac{e^{-tu}}{e^{-2\sqrt{\pi}i\sqrt{u}}}(-du)\right)
    \end{align}
    [/tex]
     
    Last edited: Nov 20, 2005
  26. Nov 22, 2005 #25

    saltydog

    User Avatar
    Science Advisor
    Homework Helper

    The third integral of the inversion is zero by inspection. More interestingly, is to show the following relation:

    [tex]
    -\frac{1}{2\pi i}\left(\int_0^{\infty}\frac{e^{-tu}}{e^{2\sqrt{\pi}i\sqrt{u}}}du-\int_{0}^{\infty}\frac{e^{-tu}}{e^{-2\sqrt{\pi}i\sqrt{u}}}du\right)=\frac{e^{-\pi/t}}{t^{3/2}}[/tex]

    Completing the square of both integrals and using the same substitutions as above:

    [tex]\sigma^2=u;\quad v=\sqrt{t}\sigma[/tex]

    I obtain:

    [tex]\int_0^{\infty}\frac{e^{-tu}}{e^{2\sqrt{\pi}i\sqrt{u}}}du=\frac{2 e^{-\pi/t}}t}\int_0^{\infty}e^{-(v+i\sqrt{\pi/t})^2}vdv[/tex]

    [tex]\int_{0}^{\infty}\frac{e^{-tu}}{e^{-2\sqrt{\pi}i\sqrt{u}}}du=\frac{2 e^{-\pi/t}}t}\int_0^{\infty}e^{-(v-i\sqrt{\pi/t})^2}vdv[/tex]

    Evaluating both of these integrals:

    [tex]\frac{2 e^{-\pi/t}}{t}\int_0^{\infty}e^{-(v+i\sqrt{\pi/t})^2}vdv=\frac{e^{-\pi/t}}{t}\left(-\frac{\pi i}{t^{1/2}}+e^{\pi/t}-\frac{\pi}{\sqrt{t}}\text{Erfi}(\sqrt{\pi/t})\right)[/tex]

    [tex]\frac{2 e^{-\pi/t}}{t}\int_0^{\infty}e^{-(v-i\sqrt{\pi/t})^2}vdv=\frac{e^{-\pi/t}}{t}\left(\frac{\pi i}{t^{1/2}}+e^{\pi/t}-\frac{\pi}{\sqrt{t}}\text{Erfi}(\sqrt{\pi/t})\right)[/tex]

    Finally, substituting these into the expression for the inverse transform:

    [tex]
    \begin{align*}
    \frac{1}{2\pi i}\int_{\gamma-i\infty}^{\gamma+i\infty}\frac{e^{zt}}{e^{2\sqrt{\pi z}}}dz &=
    -\frac{1}{2 \pi i}\left\{\frac{e^{-\pi/t}}{t}\left(-\frac{\pi i}{t^{1/2}}+e^{\pi/t}-\frac{\pi}{\sqrt{t}}\text{Erfi}(\sqrt{\pi/t})\right)\right \\ \\

    &-\frac{e^{-\pi/t}}{t}\left(\left \frac{\pi i}{t^{1/2}}+e^{\pi/t}-\frac{\pi}{\sqrt{t}}\text{Erfi}(\sqrt{\pi/t})\right)\right\} \\ \\
    &=\frac{e^{-\pi/t}}{t^{3/2}}
    \end{align}
    [/tex]

    Ok, I'm done.:smile:
     
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook