无法在 Maxima 中绘制 ODE 的解

问题描述

一天中的美好时光！

代码如下：

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,1*%pi,[ytics,false],/*[yx_ratio,0.6],*/

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

这里是错误：

integrate：变量不能是数字；发现：-12.56637061435917

出了什么问题？ 谢谢。

解决方法

看起来ode2不知道如何彻底解决问题，所以结果包含一个积分：

(%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)

我也用 contrib_ode 试过：

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

所以contrib_ode也没有完全解决。

然而，ode2 返回的解决方案（与 contrib_ode 相同）似乎是一个有效的解决方案。我将发布一个单独的答案，描述如何以数字方式评估它以进行绘图。

,

这是一种绘制解决方案 sol 的方法，该解决方案由 ode2 和 ic2 找到，如您所示。首先将 integrate 名词替换为对 quad_qags（一个数值求积函数）的调用。我将引入一个虚构的变量名称（所谓的 gensym）以避免与变量 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,epsrel = 1.0E-8,epsabs = 0.0,limit = 200)[
                 1]
                +quad_qags(%e^cos(t88336),t88336,limit = 200)[
                 1]+1)

现在我将使用该结果定义一个函数 foo1。我会列出一个数值列表，看看它是否正常工作。

(%i60) foo1(t) := ''%;
(%o60) foo1(t):=-%e
            ^((-t)-quad_qags(%e^cos(t88373),limit = 200)[
                   1]
                  +quad_qags(%e^cos(t88336),limit = 200)[
                   1]+1)
(%i61) foo1(0.5);
(%o61) -1.648721270700128
(%i62) makelist (foo1(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]

您觉得 %o62 合适吗？我会假设它没问题。接下来，我将定义一个函数 foo，当参数是数字时，它调用之前定义的 foo1，否则它只返回 0。这是 plot2d 中错误的解决方法 ，它错误地确定 foo1 不是单独的 t 函数。通常不需要该解决方法，但在这种情况下需要它。

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

好的，现在可以绘制函数 foo！

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

绘制大约需要 30 秒 - 调用 quad_qags 的成本相对较高。