2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки




На страницу 1, 2  След.
 
 Помогите в Maple и Mathematica решить систему дифуров
Сообщение20.11.2008, 13:42 
есть система:
$ 
   \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 не знаю.

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

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

 
 
 
 
Сообщение21.11.2008, 00:32 
...

 
 
 
 
Сообщение21.11.2008, 22:39 
Аватара пользователя
Код вполне боеспособен, если убрать ограничение максимального количества шагов.
Код:
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 шаг в одну тысячную, а прошагать надо до тысячи, миллион решений системы уравнений... даже не знаю сколько она это делать будет. Тут лучше уж сразу разностную схему наваять да и решать алгебраические уравнения, а не дифференциальные. Хоть мэйплом, хоть математикой, хоть чем угодно.

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

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

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

 
 
 
 
Сообщение23.11.2008, 17:28 
Код в мат лабе:
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 секунд:

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

 
 
 
 Maple.Digits
Сообщение23.11.2008, 18:51 
bvo2002 писал(а):
Может кто знает хотябы как точность решения понизить
До версии 10, количество знаков определяется значением, присвоенным переменной Digits. Более позние версии не использовал, не знаю.

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

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

 
 
 
 
Сообщение23.11.2008, 22:18 
Аватара пользователя
 !  bvo2002, для оформления текстов программ нужно использовать тэг [ code ] ... [ /code ].

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

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

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

 
 
 
 Maple
Сообщение24.11.2008, 10:50 
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);
В результате численного решения получается «накручивающаяся на окружность» траектория. Чтобы это лучше было видно использован «анимационный» вывод траектории.

 
 
 
 
Сообщение27.11.2008, 20:57 
Господа я поправил вопрос темы! Может кто теперь сможет посоветовать что или помоч!

 
 
 
 
Сообщение27.11.2008, 22:47 
Аватара пользователя
Вот надо было сразу постановку задачи писать :wink:
Математика без особых запаров решит вам такую систему (на все про все несколько секунд)
Код:
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 не особо получилось, получилась на самом деле какая-то окружность, хотя я не сильно заморачивался с тем чтобы привести все координатные оси в такой вид, чтобы они были читабельны.

 
 
 
 
Сообщение29.11.2008, 00:14 
Leierkastenmann
А если мне нужно решать методом Рунги - Кутта и погрешности должны быть определенные можно это задать?

 
 
 [ Сообщений: 20 ]  На страницу 1, 2  След.


Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group