Finding the derivative of a numeric solution to an ode in maple

Question

I'm looking at functions of the form

$w''''(t)+aw''(t)+f(w(t)) = 0$

In maple I used the following which will return numeric solutions to w,w',w'',w''':

sol := dsolve({de, ics}, type = numeric, output = listprocedure}

What I would like to be able to do is plot $w^{(4)}$ or higher order derivatives. I tried adding numeric solutions I recovered for w'' and w to return $w^{(4)}$, but I'm not quite sure how to add numeric solutions. Thanks for your help in advanced.

Answer 1

You could augment the system of DEs (including adding more ICs), so that new dummy names got associated with the higher derivatives.

That way, dsolve/numeric itself does most of the heavy numerical lifting. Probably better than trying numerical derivatives from the smaller, original result (using fdiff or evalf@D) since this way dsolve numeric can tie in any requested accuracy to however hard it might have to work to get those derivatives.

For example,

restart:

# This is your original system.
de := diff(w(t),t$4)+a*diff(w(t),t$2)+sin(w(t))=0;
ics := (D@@3)(w)(0)=1.1, (D@@2)(w)(0)=2.3, (D)(w)(0)=1.7, w(0)=Pi/2;

# This is the highest derivative you want, say.
highestdereq := diff(w(t),t$6)=d6w(t);

# Either assign to the parameters, of use the parameters=... option
# of dsolve/numeric. That's a whole other topic, so we'll keep it simple.
a:=3.2:

# An IC for the highest deriv. in the original system, which falls
# out of the original system and ICs.
eval(convert(eval(de,t=0),D),[ics]);
icw4 := isolate( %, (D@@4)(w)(0) );

# Differentiate one of the original DEs (you had just one) which
# contains the previous highest deriv. dealt with. The find its
# IC too. (This step may be repeated. Automate to get fancy.)
nextde := diff(de, t);
eval(convert(eval(nextde,t=0),D),[ics]);
icw5 := isolate( %, (D@@5)(w)(0) );

# Now invoke dsolve/numeric on the augmented system.
sol := dsolve({de, ics, highestdereq, icw4, icw5},
          type = numeric, output = listprocedure);

# If you want the procs separately... (could be done in a loop...)
W[6] := eval(d6w(t), sol):
W[5] := eval(diff(w(t),t$5), sol):
    W[4] := eval(diff(w(t),t$4), sol):
W[3] := eval(diff(w(t),t$3), sol):
    W[2] := eval(diff(w(t),t$2), sol):
W[1] := eval(diff(w(t),t$1), sol):
W[0] := eval(w(t), sol):

plot([W[0],W[1],W[2],W[3],W[4],W[5],W[6]], 0..1,
     legend=[w(t),seq('diff'(w(t),seq(t,j=1..i)),i=1..6)]);