Maple. Dsolve and functions

Question

I have a differential equation that my program solves:

p := dsolve({ic, sys}, numeric, method = rosenbrock);

After solving we have the following:

print(p(10));
      [t = 10., alpha(t) = HFloat(0.031724302221312055), beta(t) = HFloat(0.00223975915581258)]

I need use this alpha(t) and beta(t) as follows:

a := t->exp( int(alpha(t)),x=0..t) );
b := t->exp( int(beta(t)),x=0..t) )

And draw a plot:

odeplot(p, [[t, a(t)], [t, b(t)]], 0 .. 20, thickness = 2, numpoints = 500, color = [red, blue])

The first thing that occurred to do so:

p := dsolve({sys, ic},  numeric, method=rosenbrock); 
alpha := t->rhs(p(t)[2] );
beta := t->rhs(p(t)[3;
a := t->exp( int(alphat)),x=0..t) );
b := t->exp( int(betat)),x=0..t) );
odeplot(p, [[t, a(t)], [t, b(t)]], 0 .. 20, thickness = 2, numpoints = 500, color = [red, blue])

But the code does not work, and Yes, obviously, should act differently.

Answer 1

First, let's try to get your approach to work, with a few syntax and usage changes.

You didn't supply the example's details, so I make up a system of differential equations sys and initial conditions ic so that your computational commands can be performed.

restart:

sys := diff(alpha(t),t) = -1/200*beta(t),
       diff(beta(t),t) = -1/200*alpha(t) - 1/Pi^2:
ic := alpha(0)=0, beta(0)=-1:

p := dsolve({ic, sys}, numeric, method = rosenbrock, output=listprocedure):

alphat := eval(alpha(t),p):
betat := eval(beta(t),p):

a := unapply( exp( Int(alphat , 0..t) ), t, numeric):
b := unapply( exp( Int(betat , 0..t) ), t, numeric):

evalf(a(20.0)), evalf(b(20.0));

                                              -18
                   5.347592595, 3.102016550 10   

st := time():
P := plots:-odeplot(p, [[t, a(t)], [t, b(t)]], 0 .. 20,
               thickness = 2, numpoints = 50, color = [red, blue]):
( time() - st )*`seconds`;

                           16.770 seconds

P;

I used output=listprocedure so that I could assign the right procedures from solution p to alphat and betat. Those are more efficient to call many times, as opposed to your original which formed a sequence of values for each numeric t value and then had to pick off a certain operand. It's also more robust since its not sensitive to positions (which could change due to a new lexicographic ordering if you altered the names of your variables).

The above took about 16 seconds on an Intel i7. That's not fast. One reason is that the numeric integrals are being computed to higher accuracy than is necessary for the plotting. So let's restart (to ensure fair timing) and recompute with relaxed tolerances on the numeric integrations done for a and b.

restart:
sys := diff(alpha(t),t) = -1/200*beta(t),
       diff(beta(t),t) = -1/200*alpha(t) - 1/Pi^2:
ic := alpha(0)=0, beta(0)=-1:
p := dsolve({ic, sys}, numeric, method = rosenbrock, output=listprocedure):
alphat := eval(alpha(t),p):
betat := eval(beta(t),p):

a := unapply( exp( Int(alphat , 0..t, epsilon=1e-5) ), t, numeric):
b := unapply( exp( Int(betat , 0..t, epsilon=1e-5) ), t, numeric):

evalf(a(20.0)), evalf(b(20.0));

                                              -18
                   5.347592681, 3.102018090 10   

st := time():
P := plots:-odeplot(p, [[t, a(t)], [t, b(t)]], 0 .. 20,
               thickness = 2, numpoints = 50, color = [red, blue]):
( time() - st )*`seconds`;

                            0.921 seconds

You can check that this gives a plot that appears the same.

Now lets augment the example so that the numeric integrals are computed by dsolve,numeric itself. We can do that by using Calculus. In this way we are leaving it to the numeric ode solver to do its own error estimation, stepsize control, stiffness or singularity detection, etc.

Note that the integrands alphat and betat are not functions for which we have explicit functions -- their accuracy is intimately tied to the numerical ode solving. This is not quite the same as simplistically using a numerical ode routine to replace a numeric quadrature routine for a problem with an integrand which we expect to be computed directly to any desired accuracy (including on either side of any singularity).

restart:

sys := diff(alpha(t),t) = -1/200*beta(t),
       diff(beta(t),t) = -1/200*alpha(t) - 1/Pi^2,
       diff(g(t),t) = alpha(t), diff(h(t),t) = beta(t):
ic := alpha(0)=0, beta(0)=-1,
      g(0)=0, h(0)=0:

p := dsolve({ic, sys}, numeric, method = rosenbrock, output=listprocedure):

alphat := eval(alpha(t),p):
betat := eval(beta(t),p):
gt := eval(g(t),p):
ht := eval(h(t),p):

exp(gt(20.0)), exp(ht(20.0));

                                                   -18
              5.34759070530497, 3.10201330730572 10   

st := time():
P := plots:-odeplot(p, [[t, exp(g(t))], [t, exp(h(t))]], 0 .. 20,
               thickness = 2, numpoints = 50, color = [red, blue]):
( time() - st )*`seconds`;

                            0.031 seconds
P;