numericodemaplenorm

How to make a change in a numeric solution to an ODE after it has been solved?

I have two Second Order ODE's of which I want to calculate the error between them using (with Maple)

|| B(t) - A(w*t) || = sqrt(int(B(t) - A(w*t) )^2), t = 0..T)

Where A(t) is the solution to the system without any time transformation on the input and B(t) is the solution to the system with a time transformation on the input ( the time transformation is wt (w is usually small) and T is a number, not a variable. I would change T to study the error with different time spans).

For example (to help explain my question):

The original ODE is: 

diff(g(t), t, t) = -(diff(g(t), t))-sin(g(t))*cos(g(t))+sin(t)

Let A(t) be the NUMERIC solution to the original ODE (because maple cannot solve it symbolically).

Now, the ODE WITH a transformation of time in the input: 

diff(y(t), t, t) = -(diff(y(t), t))-sin(y(t))*cos(y(t))+sin(w*t)

Let B(t) be the NUMERIC solution to this ODE (because maple cannot solve it symbolically).

My question is: Is there any way to solve the error for different values of w? For this, I would need to change my numeric solution of A(t) to A(wt) after having solved for A(t) numerically. My ultimate goal is to plot the error v.s. the frequency, which is w. When w is 1, there should be no error because there is no change in the system.

I am still relatively new to coding. What I have done, since Maple cannot solve it symbolically, is solve each one numerically (with a specific w. however, I would want to do it for w in the range [0..1.5]). Then I plotted them on the same figure. However, this gives me the numeric value of A(t) NOT A(wt). And, I do not know how to subtract them. 

sol1 := dsolve([diff(g(t), t, t) = -(diff(g(t), t))- sin(g(t))*cos(g(t))+sin(t), g(0) = 0, (D(g))(0) = 0], numeric);

sol2 := dsolve([diff(y(t), t, t) = -(diff(y(t), t))-sin(y(t))*cos(y(t))+sin(.5*t), y(0) = 0, (D(y))(0) = 0], numeric);

S1 := plots[odeplot](sol1, 0 .. 10, color = red);
S2 := plots[odeplot](sol2, 0 .. 10, color = blue);
display(S1, S2);

However, this is not helpful because it is only plotting A(t), NOT A(wt). Likewise, it only plots it, it does not show me the error between them.

I am expecting the error to approach 0 as the frequency, w, approaches 0. I do expect the error to increase when w is in between 0 and 1.

Solution

  • There are a variety of ways to do this. Some are more efficient that others. Some are more convenient. I'll show a few.

    This first basic way involves two calls to dsolve, similar to your original code, but with the extra output=listprocedure option. This makes dsolve return a list of scalar-valued procedures rather than a single procedure (which would return a list of values). This allows us to pick off the individual procedures to y(t) and g(t) and utilize them separately.

    restart;
sol1 := dsolve([diff(g(t), t, t)
                = -(diff(g(t), t))- sin(g(t))*cos(g(t))+sin(t),
                g(0) = 0, (D(g))(0) = 0], numeric,
                output=listprocedure):
sol2 := dsolve([diff(y(t), t, t)
                = -(diff(y(t), t))-sin(y(t))*cos(y(t))+sin(.5*t),
                y(0) = 0, (D(y))(0) = 0], numeric,
                output=listprocedure):

    You could still use plots:-odeplot here if you wish.

    S1 := plots:-odeplot(sol1, 0 .. 20, color = red, numpoints=200):
S2 := plots:-odeplot(sol2, 0 .. 20, color = blue, numpoints=200):
plots:-display(S1, S2, size=[400,150]);

    enter image description here

    But you could also extract the individual procedures, plot them, or plot their difference.

    G := eval(g(t),sol1):
Y := eval(y(t),sol2):

plot([G,Y], 0..20, color=[red,blue], size=[400,150]);

    enter image description here

    The difference between them is now easier to plot.

    plot(G-Y, 0..20, color=black, size=[400,150]);

    enter image description here

    We can compute a norm (your integral),

    sqrt(evalf(Int( t -> ( G(t) - Y(t) )^2, 0..20, epsilon=1e-5)));
         8.74648880831735

    But now let's make it more convenient to treat w as a parameter which we adjust on the fly. (We don't want to incur the cost or inconvenience of calling dsolve for each value of w.)

    solgen := dsolve([diff(y(t), t, t)
                 = -(diff(y(t), t))-sin(y(t))*cos(y(t))+sin(w*t),
                 y(0) = 0, (D(y))(0) = 0], numeric,
                 parameters = [w],
                 output=listprocedure):

 Ygen := eval(y(t),solgen):

 # set the parameter value
 Ygen(parameters=[w=0.5]);
                       [w = 0.5]

 # we could query the parameter value
 Ygen(parameters);
                       [w = 0.5]

 # Now call it at a particular value of t
 Ygen(3.2);
                   0.864744459594622

    Now we construct a procedure of both w and t. When called without numeric arguments it returns unevaluated. When called with numeric arguments it checks whether the w value matches the current stored parameter value and if different it sets the stored value. Then it calls the computional procedure Ygen at the given t value.

    Yw := proc(t,w)
  if not [t,w]::list(numeric) then
    return 'procname'(args);
  end if;
  if w - eval(':-w',Ygen(parameters)) <> 0.0 then
    Ygen(':-parameters'=[':-w'=w]);
  end if;
  Ygen(t);
end proc:

    This can produce the same plots as before.

    plots:-display(
  plot(Yw(t,0.5), t=0..20, color=red),
  plot(Yw(t,1.0), t=0..20, color=blue),
  size=[400,150]
);

    enter image description here

    # somewhat inefficient since for each t value it
# switches the value of the w parameter.
plot(Yw(t,1.0)-Yw(t,0.5), t=0..20, size=[400,150]);

    enter image description here

    We could also do (joined lines) in point-plots, which is more efficient since all the t values are computed in sequence for any given w value. (Ie, it doesn't bounce around between w values.)

    a,b := 0,20:
xpts := Vector(100,(i)->a+(b-a)*(i-1)/(100-1),datatype=float[8]):
plots:-display(
  plot(<xpts | map(t->Yw(t,0.5), xpts)>, color=red),
  plot(<xpts | map(t->Yw(t,1.0), xpts)>, color=blue),
  size=[400,150]
);

    enter image description here

    # more efficient since all the Yw values for w=1.0 are
# computed together, and all the Yw values for w=0.5 are
# computed together.
plot(<xpts | map(t->Yw(t,1.0), xpts) - map(t->Yw(t,0.5), xpts)>,
     size=[400,150]);

    enter image description here

    evalf([seq( ['w'=w, sqrt(Int( unapply( (Yw(t,1.0)
                                         - Yw(t,w))^2, t),
                              0..20, epsilon=1e-1))],
             w=0..1.0, 0.1 )]);

    [[w = 0., 2.97123678486962], [w = 0.1, 20.3172660599286], 

     [w = 0.2, 21.8005838723429], [w = 0.3, 13.0097728518328], 

     [w = 0.4, 9.28961600039575], [w = 0.5, 8.74643983270251], 

     [w = 0.6, 6.27082651159143], [w = 0.7, 5.38965679479886], 

     [w = 0.8, 5.21680809275065], [w = 0.9, 3.19786559349464], 

     [w = 1.0, 0.]]


plot(w->sqrt(Int( (Yw(t,1.0) - Yw(t,w))^2, t=0..20,
                  epsilon=1e-3, method=_d01ajc )),
     0..1, size=[400,150]);

    enter image description here

    plot3d( Yw(t,w), w=0..1.0, t=0..50, grid=[179,179] );

    enter image description here

    # For 3D plotting it's also faster to compute all t-values
# for each given w-value, rather than the other way round.
CodeTools:-Usage(
  plot3d( Yw(t,w), w=0..1.0, t=0..50, grid=[49,49] )
):

   memory used=37.51MiB, alloc change=0 bytes, cpu time=504.00ms,
   real time=506.00ms, gc time=127.31ms

CodeTools:-Usage(
  plot3d( Yw(t,w), t=0..50, w=0..1.0, grid=[49,49] ) ):

   memory used=124.13MiB, alloc change=0 bytes, cpu time=2.66s,
   real time=2.66s, gc time=285.47ms

    [edited] The plot above of that norm is not fast. Below I make three adjustments for better performance: 1) Use a dedicated procedure Yw1 for w=1.0 so that Yw will not be called to set the parameter w twice for each evaluation of the integrand (in the norm). 2) Use option remember on that procedure Yw1. 3) Use option compile=true in the two dsolve calls.

    I also correct the formula in the norm to call Yw1(w*t), to match the original formulation in question.

    restart;
solw1 := dsolve([diff(y(t), t, t)
             = -(diff(y(t), t))-sin(y(t))*cos(y(t))+sin(t),
             y(0) = 0, (D(y))(0) = 0], numeric,
             compile=true,
             output=listprocedure):
Yw1raw := eval(y(t),solw1):
Yw1 := proc(t)
  option remember;
  if not t::numeric then return 'procname'(args); end if;
  []; # evalhf error that plot will catch
  Yw1raw(t);
end proc:
solgen := dsolve([diff(y(t), t, t)
             = -(diff(y(t), t))-sin(y(t))*cos(y(t))+sin(w*t),
             y(0) = 0, (D(y))(0) = 0], numeric,
             parameters = [w],
             compile=true,
             output=listprocedure):
Ygen := eval(y(t),solgen):
Yw := proc(t,w)
  if not [t,w]::list(numeric) then
    return 'procname'(args);
  end if;
  if w - eval(':-w',Ygen(parameters)) <> 0.0 then
    Ygen(':-parameters'=[':-w'=w]);
  end if;
  Ygen(t);
end proc:

CodeTools:-Usage(
  plot(w->sqrt(Int( (Yw1(w*t) - Yw(t,w))^2, t=0..20,
                    epsilon=1e-3, method=_d01ajc )),
       0..1, size=[400,150])
);

    memory used=0.76GiB, alloc change=205.00MiB,
    cpu time=8.22s, real time=7.78s, gc time=761.80ms

    enter image description here