Блог

Хорошего времени суток!

Вот код:

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

а вот и ошибки:

интеграция: переменная не должна быть числом; найдено: -12.56637061435917

Что пошло не так? Спасибо.

Вот способ построить решение sol , которое было найдено ode2 вами и ic2 как вы показали. Сначала замените integrate существительные вызовами на quad_qags числовую квадратурную функцию. Я введу придуманное имя переменной (так называемый генсим), чтобы избежать путаницы с переменной 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 относительно дорог.

похоже 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 ), кажется допустимым решением. Я опубликую отдельный ответ, описывающий, как оценить его численно для построения графика.