无法在Maxima中绘制ODE的解决方案

Question

一天中的好时光！

这是代码：

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]
);

这是错误：

integrate：变量不能是数字；找到：-12.56637061435917

出了什么问题？谢谢

Answer 1

这里有一种绘制解决方案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,-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)

现在我将用这个结果定义一个函数foo1。我会列出一个数值列表，看看它是否正确。

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

%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, 4*%pi], [y, -5, 5]);
plot2d: some values were clipped.
(%o64) false

这需要大约30秒的时间来绘制——调用quad_qags相对昂贵。

Answer 2

看起来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, x, t);
(%o13) [x = %c*%e^('integrate(%e^cos(t),t)-t)]

因此contrib_ode也没有完全解决这个问题。

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