Tbessel

{  TBESSEL.PDE  }

{

This example shows the integration of Bessel's Equation as a test of the

time integration capabilities of FlexPDE.

Bessel's Equation for order zero can be written as

t^2*dtt(w) + t*dt(w) + t^2*w = 0

Dividing by t^2 and avoiding the pole at t=0, we can write

dtt(w) + dt(w)/t + w = 0

FlexPDE cannot directly integrate second order equations, so we define an

auxiliary variable v=dt(w) and write a coupled pair of equations

dt(v) + v/t + w = 0

dt(w) = v

We use a dummy spatial grid of two cells and solve the equation at each node.

You can try varying the value given for ERRLIM to see how it behaves.

}

title "Integration of Bessel's Equation"

select

    ngrid=1

    errlim=1e-5  { increase accuracy to prevent accumulation of errors }

Variables

    v (threshold=0.1)

    w (threshold=0.1)

definitions

    L = sqrt(2)

    t0 = 0.001    { Start integration at t=0.001 }

Initial values    { Initialize to known values at t=t0 }

   w = 1-2.25*(t0/3)^2

   v = -0.5*t0 + 0.5625*t0*(t0/3)^2

equations

    v:  dt(v) +v/t + w = 0

    w:  dt(w) =  v

boundaries

    Region 1

       start(-L,-L) line to (L,-L) to (L,L) to (-L,L) to close

time 0.001 to 4*pi    { Exclude t=0 }

plots

   for t=0.01 by 0.01 to 0.1 by 0.1 to 1 by 1 to endtime

  history(w,bessj(0,t)) at (0,0)  as "W(t) and BESSJ0(t)"

  history(w-bessj(0,t)) at (0,0)  as "Absolute Error"

   history(v,-bessj(1,t)) at (0,0)  as "V(t) and dt(BESSJ0(t))"

   history(v+bessj(1,t)) at (0,0)  as "Slope Error"

  history(deltat)

end