Помогите в Maple и Mathematica решить систему дифуров

[bvo2002]

есть система:
$ \left\{ \begin{array}{l} \frac{dx}{dt}=-y(t)-\varphi(\sigma),\\ \frac{dy}{dt}= x(t)-\varphi(\sigma),\\ \frac{dz}{dt}=-z(t)-\varphi(\sigma), \end{array} \right. $
Причем $ \sigma = x(t)+z(t)$
А функция $ \varphi$ определяется как:
$ \left\{ \begin{array}{l} \varepsilon^3, если \sigma> \varepsilon;\\ 2(x(t)+z(t)),если \sigma\in(-\varepsilon,\varepsilon)\\ - \varepsilon^3,если \sigma<-\varepsilon$. \end{array} \right. $
Начальные условия:
x(0)=0, y(0)=-0.8165, z(0)=0.
Для начала:
\varepsilon=0.1
t = 0..1256.6371
Нужно посмотреть как разные мат средства справляются с такой системой.
Проверяю 3: MatLab, Maple12, Mathematica6.
Разобрался с MatLab. Код следующий:
Первый фаил с именем fQsys.m:

function dz = fQsys(t,z)
global myeps mymu myM
    dz = zeros(3,1); %a column vector
    phi = 0;
    if abs(z(1,1)+z(3,1)) < myeps
      phi = mymu*(z(1,1)+z(3,1));
    end
    if (z(1,1)+z(3,1)) > myeps
      phi = myM*myeps^3;
    end
    if (z(1,1)+z(3,1)) < -myeps
      phi = -myM*myeps^3;
    end
    dz(1,1) = -z(2,1)-phi;
    dz(2,1) = z(1,1)-phi;
    dz(3,1) = -z(3,1)-phi;

Второй фаил:

close all
clear all
global myeps mymu myM

myeps = 0.1; mymu = 2; myM = 1;
Dmyeps = 0.1;
T_end = 400*pi;

acc = 1e-3; 
RelTol = acc; AbsTol = acc; InitialStep = acc/10; 
options = odeset('RelTol',RelTol,'AbsTol',AbsTol,'InitialStep',InitialStep,'NormControl','on'); 

[T,z] = ode45(@fQsys,[0 T_end], [0, -(2/3)^(1/2),0], options); 
 
subplot(1,2,1);
plot(z(1:length(z),1),z(1:length(z),2));axis square; grid on;
subplot(1,2,2);
plot(T, z(1:length(z),1) + z(1:length(z),3));axis square; grid on;

Все строит все хорошо!

С Maple менее удачно. Вот код:

restart:
Digits:=4:
myeps:=0.1:
mymu:=2:
myM:=1:
with(DEtools):
F2 := proc(N,t,X,XP)                   # First order system
global myeps,mymu,myM;
local phi:=0; 
    if (abs(X[1]+X[3]) < myeps) then phi := mymu*(X[1]+X[3]);
    else if (X[1]+X[3] > myeps) then phi := myM*myeps^3;
    else if (X[1]+X[3] < -myeps) then phi := -myM*myeps^3;
         end if;
     end if;
    end if;
   XP[1] := -X[2]-phi;
   XP[2] :=  X[1]-phi;
   XP[3] := -X[3]-phi; 
end proc:
DEplot(F2,[x(t),y(t)],t=0..1256.6371,number=3,[[x(0)=0,y(0)=-0.8165,z(0)=0]],stepsize=0.1,numpoints=10500,obsrange=true,linecolor=black,thickness=1,method=rkf45,abserr=0.01,relerr=0.0001);

График конечно строит, но только зависивость x и y как матлаб не может, и с точностью плохо и зависимость от t не понятно как!!!!.
Попробовал поиграть с точностью и ничего, может еще какиенибудь параметры есть????? Помогите разобраться, в хелпе не нашел(((

С mathematica вообще беда не смог пока построить, может кто знает как???? Хотябы что то примерное???

Везде используется метод Рунге-Кутта!!! Как в mathematica не знаю.

[bvo2002]

Может кто знает хотябы как точность решения понизить, может из-за того что meple работает с большими числами поэтому так долго.
То есть чтоб вместо 0.4566774234434655655242
получать 0.45668

[Leierkastenmann]

Код математики привели бы, чтобы можно было понять что у вас не получилось и почему. Я мэйпле не спец, поэтому мне слету так не понять по коду что вообще за задача у вас такая.

[bvo2002]

...

[Leierkastenmann]

Код вполне боеспособен, если убрать ограничение максимального количества шагов.

NDSolve[{x'[t] == -y[t], 
y'[t] == x[t], 
z'[t] == -z[t], x[0] == 0, y[0] == -0.8165, 
z[0] == 0 }, {x, y, z}, {t, 0, 1256}, [color=red]MaxSteps -> 3000[/color]]; 

ParametricPlot3D[Evaluate[{x[t], y[t], z[t]} /. %], 
{t, 0, 1256}, PlotPoints -> 1000]; 

А вот если вы хотите пошагово решать систему уравнений, то в математике это можно сделать, например, вот так

soli = {};
phi = 0;
myeps = 0.1;
mymu = 2;
myM = 1;
t0 = 0.1;
x0 = 0;
y0 = -0.8165;
z0 = 0;
For[i = 0, i <= 1296/t0,
 {
  sol = NDSolve[{x'[t] == -y[t], y'[t] == x[t], z'[t] == -z[t], 
     x[i*t0] == x0, y[i*t0] == y0, z[i*t0] == z0}, {x, y, z}, {t, 
     i*t0, (i + 1)*t0}],
  phi = Which[
    First[x[(i + 1)*t0] /. sol] + 
      First[z[(i + 1)*t0] /. sol] < -myeps, -myM*myeps^3,
    First[x[(i + 1)*t0] /. sol] + First[z[(i + 1)*t0] /. sol] < myeps,
     mymu*(First[x[i*t0] /. sol] + First[z[i*t0] /. sol]),
    First[x[(i + 1)*t0] /. sol] + First[z[(i + 1)*t0] /. sol] > myeps,
     myM*myeps^3
    ],
  x0 = First[x[(i + 1)*t0] /. sol],
  y0 = First[y[(i + 1)*t0] /. sol],
  z0 = First[z[(i + 1)*t0] /. sol],
  soli = Append[soli, {x0, y0, z0}]
  }; i++
 ]

Хотя конечно думать она будет долго, если у вас по t шаг в одну тысячную, а прошагать надо до тысячи, миллион решений системы уравнений... даже не знаю сколько она это делать будет. Тут лучше уж сразу разностную схему наваять да и решать алгебраические уравнения, а не дифференциальные. Хоть мэйплом, хоть математикой, хоть чем угодно.

А еще позвольте полюбопытствовать - чем же матлабовское решение не устроило, если вы говорите, что матлаб сделал это быстро?

[bvo2002]

Просто мне нужно было сравнить решения этих 3-х покетов, с матлабом быстро разобрался там все понял, а вот на maple уже 2 недели пытался найти решение, естественно читал и help и инет шерстил. В итоге ничего кроме как делать это циклом не смог придумать. Mathmatica до этого вообще не разу в глаза не видел думал может этот покет может как матлаб, но видемо тоже придется циклом делать. Спасибо большое за помощь.

[Leierkastenmann]

Любопытство продолжает одолевать меня А можно взглянуть на матлабовский код? И узнать чистую постановку задачи, без привязки к каким-либо программам?

[bvo2002]

Код в мат лабе:
function dz = fQsys(t,z)
global myeps mymu myM

dz = zeros(3,1); %a column vector
phi = 0;

if abs(z(1,1)+z(3,1)) < myeps
phi = mymu*(z(1,1)+z(3,1));
end
if (z(1,1)+z(3,1)) > myeps
phi = myM*myeps^3;
end
if (z(1,1)+z(3,1)) < -myeps
phi = -myM*myeps^3;
end

dz(1,1) = -z(2,1)-phi;
dz(2,1) = z(1,1)-phi;
dz(3,1) = -z(3,1)-phi;

myeps = 0.1; mymu = 2; myM = 1;
Dmyeps = 0.1;
T_end = 400*pi;

acc = 1e-3;
RelTol = acc; AbsTol = acc; InitialStep = acc/10;
options = odeset('RelTol',RelTol,'AbsTol',AbsTol,'InitialStep',InitialStep,'NormControl','on');

z(1:1,1:3) = 0; z(1,2) = -(2/3)^(1/2);

Nstep = 7;
for i=0:Nstep
lz = length(z(:,1));
z0 = [z(lz,1) z(lz,2) z(lz,3)]

[T,z] = ode45(@fQsys,[0 T_end], z0, options);

plot(z(1:length(z),1),z(1:length(z),2));axis square; grid on;

myeps = myeps+Dmyeps;
end;

Добавлено спустя 10 минут 12 секунд:

Сравнить нужно эти покеты на предмет работы с такой системой!

[GAA]

bvo2002 писал(а):

Может кто знает хотябы как точность решения понизить

До версии 10, количество знаков определяется значением, присвоенным переменной Digits. Более позние версии не использовал, не знаю.

Добавлено спустя 4 минуты 25 секунд:

см., также, UseHardwareFloats.

[maxal]

bvo2002, для оформления текстов программ нужно использовать тэг [ code ] ... [ /code ].

[bvo2002]

GAA
версия 9.5 но если даже писать Digits сначала. Численое решение считается все равно не с такой точностью а в больших числах.

Добавлено спустя 55 секунд:

maxal
Спасибо буду знать.

[GAA]

1. Количество знаков при выводе результата в 10-ой версии можно установить при помощи interface. Например, если мы выполним указанный ниже код, то Maple 10 выведет значение $\pi$ с двумя знаками после запятой.

> interface(displayprecision=2);
> Digitas:= 12;
> evalf(Pi);

Я не использовал версию 9.5, но её часто ругают.

2. Classic Worksheet Maple 10

> eps:= 0.1:
> phi:= (sigma) -> piecewise(sigma < -eps, -eps^3,
      (-eps < sigma) and (sigma < eps), 2*sigma,
      sigma > eps, eps^3):
> nsol:= dsolve([
        diff(x(t), t) = -y(t) - phi(x(t)+z(t)),
        diff(y(t), t) =  x(t) - phi(x(t)+z(t)),
        diff(z(t), t) = -z(t) - phi(x(t)+z(t)),
        x(0)= 0, y(0) = -0.8165, z(0) = 0], numeric):
> with(plots):
> odeplot(nsol, [x(t), y(t), z(t)], frames=20);

В результате численного решения получается «накручивающаяся на окружность» траектория. Чтобы это лучше было видно использован «анимационный» вывод траектории.

[bvo2002]

Господа я поправил вопрос темы! Может кто теперь сможет посоветовать что или помоч!

[Leierkastenmann]

Вот надо было сразу постановку задачи писать
Математика без особых запаров решит вам такую систему (на все про все несколько секунд)

sol = NDSolve[
  {
   x'[t] == -y[t] + 
     Which[x[t] + z[t] < -eps, -eps^3, x[t] + z[t] > eps, eps^3, 
      x[t] + z[t] < eps, 2*(x[t] + z[t])],
   y'[t] == 
    x[t] + Which[x[t] + z[t] < -eps, -eps^3, x[t] + z[t] > eps, eps^3,
       x[t] + z[t] < eps, 2*(x[t] + z[t])],
   z'[t] == -z[t] + 
     Which[x[t] + z[t] < -eps, -eps^3, x[t] + z[t] > eps, eps^3, 
      x[t] + z[t] < eps, 2*(x[t] + z[t])],
   x[0] == 0, y[0] == -0.8165, z[0] == 0
   },
  {x, y, z},
  {t, 0, 1256}
  ]

Главное eps только численно задайте, а потом делайте ParametricPlot3D. Правда у меня 3D не особо получилось, получилась на самом деле какая-то окружность, хотя я не сильно заморачивался с тем чтобы привести все координатные оси в такой вид, чтобы они были читабельны.

[bvo2002]

Leierkastenmann
А если мне нужно решать методом Рунги - Кутта и погрешности должны быть определенные можно это задать?