Search code examples
plotintegrationmapleexpansion

how to avoid singularity in numerical integration?


I would like to plot E2(t) in Maple18, but it blows up around t=0. How can I apply the asymptotic evaluation for this integral to avoid the singularity and solve this problem?

restart: with(plots):
with(Student[NumericalAnalysis]):
g1:=(x,t)->(-sqrt(t)/(2*sqrt(Pi*r^3)))*(sin((r*(x-1)^(2)/(4*t))+(Pi/4))-sin((r*(x+1)^(2)/(4*t))+(Pi/4))):

g2:=(x,t)->int(g1(x,t),r=1..infinity);

g3:=(x,t) -> (diff(g2(x,t),t)):


g4:=(x,t) -> (diff(g2(x,t),x,x)):


g5:=(x,t) -> ((1/2)*(g3(x,t)^2+g4(x,t)^2)):


E2:=t->(int(g5(x,t),x=0..100)):

evalf(E2(0));
Error, (in g1) numeric exception: division by zero
evalf(E2(1));
Error, (in g3) invalid input: diff received 1, which is not valid for its 2nd argument
plot(E2(t),t=0..20);

Best regards,


Solution

  • I don't know what you mean by "avoid the singularity". Your function is singular at t=0, and there's no way to change that. Just don't try to evaluate it at 0.

    Your other problem is caused by trying to differentiate with respect to a number. What you need to do is create the procedure E2 with unapply. Then the differentiations will be done symbolically. Here's your code, essentially. I removed superfluous parentheses, the (x,t)->, and the unneeded package loading.

    restart:
    g1:= -sqrt(t/Pi/r^3)/2 * 
         (sin(r*(x-1)^2/4/t+Pi/4)-sin(r*(x+1)^2/4/t+Pi/4)):
    g2:= int(g1, r= 1..infinity):
    g3:= diff(g2,t):
    g4:= diff(g2,x$2):
    g5:= (g3^2+g4^2)/2:
    

    And here's E2:

    E2:= unapply(Int(g5, x= 0..100, epsilon= 1e-4, digits= 7), t):
    

    Note the capital I in Int. This will avoid wasting time trying symbolic integration and will go directly to numeric integration when evalf is applied. The epsilon and digits arguments will decrease the precision and increase the speed. The precision will still be adequate for a plot. The plot command is

    plot(E2, 0..20, numpoints= 50, labels= [t, ``]);
    

    This will take two to three minutes. If you increase the numpoints for a more-accurate plot, the time will increase proportionally.