Advanced Search
All Posts : Applications, Examples and Libraries
This post in reply to the Question, solve ordinary differential equations

 

While generating a 3D plot of the solution of an ODE with a parameter, I noticed that better performance could be obtained by calling the plot3d command with a procedure argument, done in a special manner.

I don't recall this being discussed before, so I'll share it. (It it has been discussed, and this is widely known, then I apologize.)

I tweaked the initial conditions of the original ODE system, to obtain a non-trivial solution. I don't think that the particular nature of the solution has a bearing on this note.

> 

restart;
> 

Digits := 15:
> 

 

 

The ODE system has two parameters. One, A, will get a fixed value. The

other, U0, will be used like an independent variable for plotting.

 

> 

 

> 

eq1:= diff(r[1](t), t, t)+(.3293064114+209.6419478*U0)*(diff(r[1](t), t))
      +569.4324330*r[1](t)-0.3123434112e-2*V(t) = -1.547206836*U0^2*q(t):
eq2:= 2.03*10^(-8)*(diff(V(t), t))+4.065040650*10^(-11)*V(t)
      +0.3123434112e-2*(diff(r[1](t), t)) = 0:
eq3 := diff(q(t), t, t)+1047.197552*U0*(q(t)^2-1)*(diff(q(t), t))
       +1.096622713*10^6*U0^2*q(t) = -2822.855019*A*(diff(r[1](t), t, t)):
> 

 

> 

ics:=r[1](0)=0,D(r[1])(0)=1e-1,V(0)=0,q(0)=0,D(q)(0)=0:
> 

 

> 

res := dsolve({eq1,eq2,eq3, ics},numeric,output=listprocedure,parameters=[A,U0]):

 

I will call the procedure returned by dsolve, for evalutions of V(t), as the

dsolve numeric solution-procedure in the discussion below.

 

> 

WV := eval(V(t), res):
> 

 

> 

WV(parameters=[A=1e0]):
> 

 

 

The goal is to produce a 3D plot of V(t) as a function of t and the parameter U0.

 

> 

 

> 

tlo,thi := 0.0, 2.0;
U0lo,U0hi := 1e-3, 2e-1;

0., 2.0

0.1e-2, .2

 

This is the grid size used for plot3d below. It is nothing special.

 

> 

(m,n) := 51,51;

51, 51

 

First, I'll demonstrate that a 3D plot can be produced quickly, by populating a
Matrix for floating-point evaluations of V(t), depending on t in the first
Matrix-dimension and on parameter U0 in the second Matrix-dimansion.

 

The surfdata command is used. This is similar to how plot3d works.

 

This  computes reasonably quickly.

 

But generating the numeric values for U0 and t , based on the i,j positions

in the Matrix, involves the kind of sequence generation formulas that are

error prone for people.

 

> 

str := time[real]():
M:=Matrix(m,n,datatype=float[8]):
for j from 1 to n do
  u0 := U0lo+(j-1)*(U0hi-U0lo)/(n-1);
  WV(parameters=[U0=u0]);
  for i from 1 to m do
    T := tlo+(i-1)*(thi-tlo)/(m-1);
    try
      M[i,j] := WV(T);
    catch:
      # mostly maxfun complaint for t above some value.
      M[i,j] := Float(undefined);
    end try;
  end do:
end do:
plots:-surfdata(M, tlo..thi, U0lo..U0hi,
                labels=["t","U0","V(t)"]);
(time[real]()-str)*'seconds';

1.686*seconds

 

So let's try it using the plot3d command directly. A 2-parameter procedure
is constructed, to pass to plot3d. It's not too complicated. This procedure
will uses one of its numeric arguments to set the ODE's U0 parameter's

value for the dsolve numeric solution-procedure, and then pass along

the other numeric argument as a t value.


It's much slower than the surfdata call above..

 

> 

VofU0 := proc(T,u0)
       WV(parameters=[U0=u0]);
       WV(T);
     end proc:
> 

str := time[real]():
plot3d(VofU0, tlo..thi, U0lo..U0hi,
       grid=[m,n], labels=["t","U0","V(t)"]);
(time[real]()-str)*'seconds';

37.502*seconds

 

One reason why the previous attempt is slow is that the plot3d command

is changing values for U0 in its outer loop, and changing values of t in its

inner loop. The consequence is that the value for U0 changes for every

single evaluation of the plotted procedure. This makes the dsolve numeric

solution-procedure work harder, by losing/discarding prior numeric
solution details.

 

The simple 3D plot below demonstrates that the plot3d command chooses
x-y pairs by letting its second supplied independent variable be the one
that changes in its outer loop. Each time the value for y changes the counter

goes up by one.

 

> 

glob:=0:
plot3d(proc(x,y) global glob; glob:=glob+1; end proc,
       0..3, 0..7, grid=[3,3],
       shading=zhue,  labels=["x","y","glob"]);

 

So now let's try and be clever and call the plot3d command with the two
independent variables reversed in position (in the call). That will make

the outer loop change t instead of the ODE parameter U0.

 

We can use the transform command to swap the two indepenent
axes in the plot, if we prefer the axes roles switched. Or we could use the
parametric calling sequence of plot3d for the same effect.

 

The problem is that this is still much slower!

 

> 

VofU0rev := proc(u0,T)
       WV(parameters=[U0=u0]);
       WV(T);
     end proc:
> 

str := time[real]():
Prev:=plot3d(VofU0rev, U0lo..U0hi, tlo..thi,
             grid=[n,m], labels=["U0","t","V(t)"]):
(time[real]()-str)*'seconds';

plots:-display(
  plottools:-transform((x,y,z)->[y,x,z])(Prev),
  labels=["t","U0","V(t)"],
  orientation=[50,70,0]);

34.306*seconds

 

There is something else to adjust, to get the quick timing while using

the plot3d command here.

 

It turns out that setting the parameter's numeric value in the
dsolve numeric solution-procedure causes the loss of previous details
of the numeric solving, even if the parameter's value is the same.

 

So calling the dsolve numeric solution-procedure to set the parameter

value must be avoided, in the case that the new value is the same as

the old value.

 

One way to do that is to have the plotted procedure first call the

dsolve numeric solution-procedure to query the current parameter

value, so as to not reset the value if it is not changed. Another way

is to use a local of an appliable module to store the running value

of the parameter, and check against that. I choose the second way.

 

And plot3d must still be called with the first independent variable-range

as denoting the ODE's parameter (instead of the ODE's independent

variable).

 

And the resulting plot is fast once more.

 

> 

VofU0module := module()
       local ModuleApply, paramloc;
       ModuleApply := proc(par,var)
         if not (par::numeric and var::numeric) then
           return 'procname'(args);
         end if;
         if paramloc <> par then
           paramloc := par;
           WV(parameters=[U0=paramloc]);
         end if;
         WV(var);
       end proc:
end module:

 

For fun, this time I'll use the parameter calling sequence to flip the

axes, instead of plots:-transform. That's just because I want t displayed

on the first axis. But for the performance gain, what matters is that it

is U0 which gets values from the first axis plotting-range.

 

> 

str := time[real]():
plot3d([y,x,VofU0module(x,y)], x=U0lo..U0hi, y=tlo..thi,
       grid=[n,m], labels=["t","U0","V(t)"]);
(time[real]()-str)*'seconds';

1.625*seconds

 

And, naturally, I could also use the parametric form to get a fast plot

with the axes roles switched.

 

> 

str := time[real]():
plot3d([x,y,VofU0module(x,y)], x=U0lo..U0hi, y=tlo..thi,
       grid=[n,m], labels=["U0","t","V(t)"]);
(time[real]()-str)*'seconds';

1.533*seconds

> 

 

Download ode_param_plot.mw
ode parameter dsolve numeric plotting

Please Wait...