Unable to plot solution of ODE in Maxima

Question

Good time of the day!

Here is the code:

eq:'diff(x,t)=(exp(cos(t))-1)*x;
ode2(eq,x,t);
sol:ic1(%,t=1,x=-1);

/*---------------------*/

plot2d(

rhs(sol), 

[t,-4*%pi,  4*%pi],
[y,-5,5],
[xtics,-4*%pi, 1*%pi, 4*%pi],

[ytics, false], 

/*[yx_ratio , 0.6], */

[legend,"Solution."],
[xlabel, "t"], [ylabel, "x(t)"],
[style, [lines,1]],
[color, blue]
);

and here is the errors:

integrate: variable must not be a number; found: -12.56637061435917

What went wrong? Thanks.

Answer 1

Here's a way to plot the solution sol which was found by ode2 and ic2 as you showed. First replace the integrate nouns with calls to quad_qags, a numerical quadrature function. I'll introduce a made-up variable name (a so-called gensym) to avoid confusion with the variable t.

(%i59) subst (nounify (integrate) = 
              lambda ([e, xx], 
                      block ([u: gensym(string(xx))], 
                             quad_qags (subst (xx = u, e), u, -4*%pi, xx)[1])),
              rhs(sol)); 
(%o59) -%e^((-t)-quad_qags(%e^cos(t88373),t88373,-4*%pi,t,
                           epsrel = 1.0E-8,epsabs = 0.0,
                           limit = 200)[
                 1]
                +quad_qags(%e^cos(t88336),t88336,-4*%pi,t,
                           epsrel = 1.0E-8,epsabs = 0.0,
                           limit = 200)[
                 1]+1)

Now I'll define a function foo1 with that result. I'll make a list of numerical values to see if it works right.

(%i60) foo1(t) := ''%;
(%o60) foo1(t):=-%e
            ^((-t)-quad_qags(%e^cos(t88373),t88373,-4*%pi,t,
                             epsrel = 1.0E-8,epsabs = 0.0,
                             limit = 200)[
                   1]
                  +quad_qags(%e^cos(t88336),t88336,-4*%pi,t,
                             epsrel = 1.0E-8,epsabs = 0.0,
                             limit = 200)[
                   1]+1)
(%i61) foo1(0.5);
(%o61) -1.648721270700128
(%i62) makelist (foo1(t), t, makelist (k, k, -10, 10));
(%o62) [-59874.14171519782,-22026.46579480672,
        -8103.083927575384,-2980.957987041728,
        -1096.633158428459,-403.4287934927351,
        -148.4131591025766,-54.59815003314424,
        -20.08553692318767,-7.38905609893065,-2.71828182845904,
        -1.0,-0.3678794411714423,-0.1353352832366127,
        -0.04978706836786394,-0.01831563888873418,
        -0.006737946999085467,-0.002478752176666358,
        -9.118819655545163E-4,-3.354626279025119E-4,
        -1.234098040866796E-4]

Does %o62 look right to you? I'll assume it is okay. Next I'll define a function foo which calls foo1 defined before when the argument is a number, otherwise it just returns 0. This is a workaround for a bug in plot2d, which incorrectly determines that foo1 is not a function of t alone. Usually that workaround isn't needed, but it is needed in this case.

(%i63) foo(t) := if numberp(t) then foo1(t) else 0;
(%o63) foo(t):=if numberp(t) then foo1(t) else 0

Okay, now the function foo can be plotted!

(%i64) plot2d (foo, [t, -4*%pi, 4*%pi], [y, -5, 5]);
plot2d: some values were clipped.
(%o64) false

That takes about 30 seconds to plot -- calling quad_qags is relatively expensive.

Answer 2

it looks like ode2 does not know how to completely solve the problem, so the result contains an integral:

(%i6) display2d: false $
(%i7) eq:'diff(x,t)=(exp(cos(t))-1)*x;
(%o7) 'diff(x,t,1) = (%e^cos(t)-1)*x
(%i8) ode2(eq,x,t);
(%o8) x = %c*%e^('integrate(%e^cos(t),t)-t)
(%i9) sol:ic1(%,t=1,x=-1);
(%o9) x = -%e^((-%at('integrate(%e^cos(t),t),t = 1))
              +'integrate(%e^cos(t),t)-t+1)

I tried it with contrib_ode also:

(%i12) load (contrib_ode);
(%o12) "/Users/dodier/tmp/maxima-code/share/contrib/diffequations/contrib_ode.mac"
(%i13) contrib_ode (eq, x, t);
(%o13) [x = %c*%e^('integrate(%e^cos(t),t)-t)]

So contrib_ode did not solve it completely either.

However the solution returned by ode2 (same for contrib_ode) appears to be a valid solution. I'll post a separate answer describing how to evaluate it numerically for plotting.