Maple Аппроксимация решений систем дифф.уравнений

Форум пользователей пакета Maple

Модератор: Admin

teolog
Сообщения: 111
Зарегистрирован: Вс янв 22, 2012 8:41 pm

Maple Аппроксимация решений систем дифф.уравнений

Сообщение teolog » Вс янв 22, 2012 8:46 pm

Здравствуйте. Нигде не могу найти хотябы намёк на решение задачи. Может вы могли бы чем-либо помочь?
Дана система дифференциальных уравнений в Maple:

Код: Выделить всё

> sys:=diff(x1(t),t)=-4*x1(t)^2, diff(x2(t),t)=x1(t)^2-x2(t);
> dsolve({sys,x1(0)=x10, x2(0)=x20},{x1(t),x2(t)});

Решением данной системы будет

Код: Выделить всё

{x1(t) = 1/(4*t+1/x10), x2(t) = (-(1/16)*exp(t)/(t+1/(4*x10))-(1/16)*exp(-1/(4*x10))*Ei(1, -t-1/(4*x10))+(1/4)*x10+(1/16)*exp(-1/(4*x10))*Ei(1, -1/(4*x10))+x20)*exp(-t)}

Ei - интегральная показательная функция, которая относится к специальным функциям, интегралы от которых не берутся. Т.е. эта функция не может быть выражена через элементарные функции. Сам Maple без труда строит график решения (исключается параметр времени, и строится двумерная кривая х2(х1) как функция одной координаты от другой):

Код: Выделить всё

> with(plots):
> sys:= diff(x1(t),t)=-4*x1(t)^2, diff(x2(t),t)=x1(t)^2-x2(t);
> fens:={x1(t),x2(t)}:
> F:=dsolve({sys, x1(0)=1, x2(0)=1}, fens, type=numeric):
> odeplot(F,[x1(t),x2(t)],0..5,numpoints=200);

Вопрос. как можно аппроксимировать аналитической функцией графическое решение? Или хотя бы как "достать" из этой кривой массив точек?

Markiyan Hirnyk
Сообщения: 1366
Зарегистрирован: Вс дек 04, 2011 11:07 pm

Ссылка

Сообщение Markiyan Hirnyk » Вс янв 22, 2012 9:01 pm


teolog
Сообщения: 111
Зарегистрирован: Вс янв 22, 2012 8:41 pm

Сообщение teolog » Вс янв 22, 2012 11:41 pm

Спасибо огромное!
а есть ли какие-то функции, чтобы сразу аппроксимирующую функцию в аналитическом виде из решения системы ДУ получить?

Markiyan Hirnyk
Сообщения: 1366
Зарегистрирован: Вс дек 04, 2011 11:07 pm

Сообщение Markiyan Hirnyk » Пн янв 23, 2012 12:02 am

teolog писал(а):Спасибо огромное!
а есть ли какие-то функции, чтобы сразу аппроксимирующую функцию в аналитическом виде из решения системы ДУ получить?
Вопрос не понял, поэтому ответить не могу. Пожалуйста приведите пример к рассматриваемой Вами системе ОДУ.

teolog
Сообщения: 111
Зарегистрирован: Вс янв 22, 2012 8:41 pm

Сообщение teolog » Пн янв 23, 2012 12:27 pm

Я Maple начал изучать недавно, и приходится обучаться походу, потому что сроки поджимают. Используя ссылку, которую Вы мне дали у меня последовательность действий следующая (хотя я наверняка делаю что-то неправильно):

Код: Выделить всё

> sys:= diff(x1(t),t)=-4*x1(t)^2, diff(x2(t),t)=x1(t)^2-x2(t);
> inicond:=x1(0)=1, x2(0)=1;
> fens:={x1(t),x2(t)};
> RES:=dsolve({sys,inicond},fens,type=numeric,output=array([seq(j/10,j=0..70)])):
> B := convert(rhs(op([1, 3, 2],RES)), array);
> writedata("d:/testfile3.txt", B);

(кстати, не могу понять что за конструкция op([1, 3, 2] и почему цифры 1 3 2 именно в таком порядке идут)
Далее открываю файл testfile3.txt, копирую оттуда данные, вставляю в Exel. вставляется в виде таблицы. Удаляю вручную лишние строки (из 70 значений, мне вполне хватит и 10-15, просто точки располагаются неравномерно, шаг не фиксированый). Транспонирую эту таблицу в строки. Эти две строки вставляю в Word. Там делаю автоматическую замену табуляции на запятые. Затем эти две строки вставляю снова в Maple, и аппроксимирую

Код: Выделить всё

> pnts:=[первая строка],[вторая строка];
> g:=interp(pnts, x); #интерполяция полиномом
> with(stats):
> h:=fit[leastsquare[[x,y], y=a1*x^3+a2*x^2+a3*x+a4]]([pnts]); #интерполяция по методу наименьших квадратов

В итоге получаю две аппроксимации решения исходного ДУ

Код: Выделить всё

g := -34814.65534*x^10+1.341932387*10^5*x^9-2.174505090*10^5*x^8+1.936324135*10^5*x^7-1.035120647*10^5*x^6+33784.43480*x^5-6394.968387*x^4+558.6577441*x^3+5.218646791*x^2-.7659385706*x+6.*10^(-10)
h := y = .8785963782*x^3-3.039856567*x^2+3.233695213*x-0.8042632574e-1

Можно ли получить такие аппроксимации сразу же? Т.е. чтобы записать точки, через которые проходит решение ДУ - потребовалось всего две короткие строчки кода. Можно ли так же, какой-либо функцией, сразу получить аппроксимированное решение ДУ? просто мне надо делать анализ очень большого числа аппроксимаций, и такая последовательность действий сильно тормозит работу, и я искал как бы её ускорить. Заранее спасибо.

Markiyan Hirnyk
Сообщения: 1366
Зарегистрирован: Вс дек 04, 2011 11:07 pm

Сообщение Markiyan Hirnyk » Пн янв 23, 2012 11:10 pm

Если я правильно понимаю Ваш вопрос, Вам после решения системы ДУ требуется явная зависимость x2(t) oт x1(t), с которой можно работать, т.е. считать, рисовать график, интегрировать, дифференцировать и. т. далее. В Вашем примере это сделать нетрудно.
> sys := diff(x1(t), t) = -4*x1(t)^2, diff(x2(t), t) = x1(t)^2-x2(t);
> G := dsolve({sys, x1(0) = 1, x2(0) = 1}, {x1(t), x2(t)});

{x1(t) = 1/(4*t+1), x2(t) = (-(1/16)*exp(t)/(t+1/4)-(1/16)*exp(-1/4)*Ei(1, -t-1/4)+(1/16)*(20*exp(1/4)+Ei(1, -1/4))/exp(1/4))*exp(-t)}

Далее, исключаем t. Поскольку с незвестными, которые являются функциями, работать затруднительно, то переобозначим их:
H := eval(G, [x1(t) = x, x2(t) = y])
{x = 1/(4*t+1), y = (-(1/16)*exp(t)/(t+1/4)-(1/16)*exp(-1/4)*Ei(1, -t-1/4)+(1/16)*(20*exp(1/4)+Ei(1, -1/4))/exp(1/4))*exp(-t)}
Если Вы желаете узнать свойства команды, например, eval, Вы набираете в рабочем поле эту команду с вопросительным знаком перед ней ?eval , наводите на нее в любом месте курсор и нажимаете клавишу Enter. Появлятся страница справки о команде с примерами и ссылками. Все изложено на английском языке. Кроме того, начальніе сведения о Maple неплохо изложены в книге Б. Манзона по адресу http://www.exponenta.ru/soft/Maple/mans ... tion/0.asp . Итак, продолжим. Находим t через x1(t) из первого уравнения системы H[1]:
solve(H[1], t)
-(1/4)*(x-1)/x
Подставляем это выражение во второе уравнение H[2] ипреобразуем выражение в функцию (см. за подробностями ?rhs ?solve ?unapply):
y := unapply(rhs(eval(H[2], t = solve(H[1], t))), x)
x -> (-(1/16)*exp(-(1/4)*(x-1)/x)/(-(1/4)*(x-1)/x+1/4)-(1/16)*exp(-1/4)*Ei(1, (1/4)*(x-1)/x-1/4)+(1/16)*(20*exp(1/4)+Ei(1, -1/4))/exp(1/4))*exp((1/4)*(x-1)/x)
Теперь можем нарисовать график функции
plot(y)
Изображение
вычислить ее значение
evalf(y(30))
-6.135101292+1.273369666*10^(-10)*I
Очень маленькую мнимую часть, возникшую в результате приближенных вычислений, можно убрать
fnormal(evalf(y(30)))
-6.135101292+0.*I
Как Вы правильно заметили, y(x)и ее первообразная не являются элементарными функциями, однако определенный интеграл можно найти численно
int(y(x), x = 1 .. 2, numeric)
.9502666328+3.392661131*10^(-14)*I
Можем также найти производную
diff(y(x), x)
(-(1/16)*(-1/(4*x)+(1/4)*(x-1)/x^2)*exp(-(1/4)*(x-1)/x)/(-(1/4)*(x-1)/x+1/4)+(1/16)*exp(-(1/4)*(x-1)/x)*(-1/(4*x)+(1/4)*(x-1)/x^2)/(-(1/4)*(x-1)/x+1/4)^2+(1/16)*exp(-1/4)*(1/(4*x)-(1/4)*(x-1)/x^2)*exp(-(1/4)*(x-1)/x+1/4)/((1/4)*(x-1)/x-1/4))*exp((1/4)*(x-1)/x)+(-(1/16)*exp(-(1/4)*(x-1)/x)/(-(1/4)*(x-1)/x+1/4)-(1/16)*exp(-1/4)*Ei(1, (1/4)*(x-1)/x-1/4)+(1/16)*(20*exp(1/4)+Ei(1, -1/4))/exp(1/4))*(1/(4*x)-(1/4)*(x-1)/x^2)*exp((1/4)*(x-1)/x)

teolog
Сообщения: 111
Зарегистрирован: Вс янв 22, 2012 8:41 pm

Сообщение teolog » Вт янв 24, 2012 1:04 am

Уважаемый Markiyan Hirnyk, не хочу Вас обидеть или показаться грубым, но и мепловскую справку, и книгу Манзона я открываю при поиске решения той или иной задачи. Также гуглю постоянно. Но конечно же, хватая кусками, системных знаний по командам и их синтаксису не приобретаю, и порой решение задачи (в смысле поиска команд и функций) существует, и весьма простое, но совершенно для меня неочевидное.
Да, Вы правы, мне нужна явная зависимость x2(t) oт x1(t). Но напишу исходную задачу, возможно Вы дадите ещё какой-либо совет для упрощения решения.
Изначально задача выглядит так. Дана система дифференциальных уравнений, описывающих динамическую систему.
> sys := diff(x1(t), t) = U, diff(x2(t), t) = x1(t)^2-x2(t);
Статическая характеристика системы представляет собой параболу
x2=x1^2
Система начинает и заканчивает движение на статической характеристике. Т.е. начальная и конечная точка принадлежат параболе.Опуская теорию оптимальных управлений скажу что я рассматриваю случай, когда система движется при управлениях
U=-1,
U=-4*x1(t)^2,
U=+1

В процессе движения от начальной точки до конечной система может быть как под действием всех трёх управлений по очереди, так и только двух (или даже одного).
Моей задачей является нахождение минимального интеграла координаты x1(t) по времени. Т.е.
int(x1(t),t=0..T)
где T не задана, а определяется лишь по достижении конечной точки.
При движении системы под действием управлений +1 и -1 система дифференциальных уравнений решается в общем виде, функции получаются аналитическими, я исключаю время, и нахожу траекторию движения системы как функцию x2(x1).
Однако пропытка решить данную систему ДУ при управлении равном U=-4*x1(t)^2 возникает неберущийся интеграл.Поэтому я попытался аппроксимировать решение системы. Тот способ который Вы предложили - очень помог мне. Однако для нахождения итоговой аппроксимированной зависимости мне приходится использовать Exel и Word, чтобы найденный Меплом массив точек снова загнать в Мепл и аппроксимировать его. Извините, но плохо знаю синтаксис (даже где это приблизительно можно искать), поэтому так ламерски длинно решаю задачу.
Изображение
Исходная задача выглядит так как на рисунке. Задавая меняющиеся начальные условия в качестве цикла, я двигаю кривую при управлении U=-4*x1(t)^2 и ищу минимальный интеграл int(x1(t),t=0..T) в пределах первого квадранта декартовой сетки. Начальная и конечная точка фиксированы, верхняя и нижняя траектории тоже. Двигается только та что между ними (средняя). Направление движения системы указано стрелками. Найти такую траекторию движения, чтобы некоторая интегральная функция была минимальной.
Для примера поясню как я решал задачу, если вместо управления U=-4*x1(t)^2 будет управление U=-x1(t). Последовательность управлений такая же: сначала U=-1, затем U=-x1(t), затем U=+1.
До этого я задачу решал так как в программе ниже. В ней x10 и x20 - начальная точка движения системы. x1t и x2t - конечная точка движения системы. x11 и x21 - точка первого переключения управления. x12 и x22 - точка второго переключения управления.

x21:=x11^2+2*x11+2+exp(x11-x10)*(-x10^2-2*x10-2+x20) -Траектория движения системы под управлением -1
x2t=x1t^2-2*x1t+2+exp(-x1t+x12p)*(-x12p^2+2*x12p-2+x22p) - Траектория движения системы под управлением +1
x22p=-x12p^2+(x11+x21/x11)*x12p - Траектория движения системы под управлением -x1.

Программа
restart;
n1:=1; n2:=1; k1:=0.9; k2:=0.81; h:=-0.05;
#Задание начальной и конечной точки для системы и шага
for x11 from n1 by h while (x11>0) do
x10:=n1; x20:=n2; x1t:=k1; x2t:=k2;
x21:=x11^2+2*x11+2+exp(x11-x10)*(-x10^2-2*x10-2+x20);
#Траектория движения системы под управлением -1
#Нахожу пересечение второй и третьей траекторий
sys1:={x2t=x1t^2-2*x1t+2+exp(-x1t+x12p)*(-x12p^2+2*x12p-2+x22p), x22p=-x12p^2+(x11+x21/x11)*x12p}:
x2:=fsolve(sys1,{x12p, x22p}, x12p=0..1);
x12:=rhs(x2[1]);
T1:=-x11+x10; T2:=-ln(x12/x11); T3:=x1t-x12;
J1:=int(-t+x10,t=0..T1)+int(x11*exp(-t), t=0..T2)+int(t+x12, t=0..T3);

f10:=n1; f20:=n2; f1t:=k1; f2t:=k2;
f11:=x11+h;
f21:=f11^2+2*f11+2+exp(f11-f10)*(-f10^2-2*f10-2+f20);
sys3:={f2t=f1t^2-2*f1t+2+exp(-f1t+f12p)*(-f12p^2+2*f12p-2+f22p), f22p=-f12p^2+(f11+f21/f11)*f12p}:
f2:=fsolve(sys3,{f12p, f22p}, f12p=0..1);
f12:=rhs(f2[1]);
t1:=-f11+f10; t2:=-ln(f12/f11); t3:=f1t-f12;
J2:=int(-t+f10,t=0..t1)+int(f11*exp(-t), t=0..t2)+int(t+f12, t=0..t3);

if (J2>J1) then
print ("x11=",x11,J1,"x12=",x12);
break;
end if;
end do;

Извините, если длинно и не очень понятно, но таким образом я искал минимум интеграла одной координаты по времени, при движении целой системы.
Но вот управление приняло более сложный вид, и таким образом приходится вместо задания цикла каждый раз вбивать новый вид второй траектории (коэффициенты при аппроксимации меняются). И вручную считать точки пересечения и интеграл. Причём ещё раз повторю, система движется в координатах x2(x1), а найти надо минимум интеграла одной координаты по времени int(x1(t),t=0..T).
Т.е. есть ли какая-либо команда/функция, чтобы из решения системы ДУ (даже сложного вида) можно было бы сразу получить аппроксимирующую функцию в виде зависимостb x2(t) oт x1(t)?
Заранее спасибо за советы.

teolog
Сообщения: 111
Зарегистрирован: Вс янв 22, 2012 8:41 pm

Сообщение teolog » Вт янв 24, 2012 1:18 am

Если ещё более обще, то может Мепл для заданного диапазона измениния x и y автоматически найти точки пересечения двух или более траекторий (траектории являются решениями систем ДУ), и посчитать вдоль этих траекторий интеграл какой-либо функции как по координатам, так и по исключённому параметру (по времени)?

Markiyan Hirnyk
Сообщения: 1366
Зарегистрирован: Вс дек 04, 2011 11:07 pm

Сообщение Markiyan Hirnyk » Чт янв 26, 2012 9:15 pm

teolog писал(а):Если ещё более обще, то может Мепл для заданного диапазона измениния x и y автоматически найти точки пересечения двух или более траекторий (траектории являются решениями систем ДУ), и посчитать вдоль этих траекторий интеграл какой-либо функции как по координатам, так и по исключённому параметру (по времени)?
Например, таким образом:
sys := diff(x1(t), t) = -4*x1(t)^2, diff(x2(t), t) = x1(t)^2-x2(t):
G := dsolve({sys, x1(0) = 1, x2(0) = 1}, {x1(t), x2(t)}):#одна интегральная кривая
evalc(Re(rhs(G[2])));
#избавляемся от мнимой части в решении

(-(1/16)*exp(t)/(t+1/4)+(1/16)*exp(-1/4)*Ei(t+1/4)+((5/4)*exp(1/4)-(1/16)*Ei(1/4))/exp(1/4))*exp(-t)

sys1 := diff(x1(s), s) = -1, diff(x2(s), s) = x1(s)^2-x2(s):
H := dsolve({sys1, x1(0) = 1, x2(0) = 1}, {x1(s), x2(s)}):#другая интегральная кривая

Digits := 20; DirectSearch:-SolveEquations([x = rhs(G[1]), y = evalc(Re(rhs(G[2]))), x = rhs(H[1]), y = rhs(H[2])], {s >= -1, -1 <= t, 0 <= x, 0 <= y, s <= .1, t <= .1, x <= 2, y <= 2});#пересечение этих интегральных кривых
[
[4.0971460009758634691 10:(-16) , Vector[column]([[3.43660448424 10^(-9)],[-1.492444323079 10^(-8)],[6.258050776 10^(-10)],[1.322020090514 10^(-8)]]), [

s = 0.00019374013908754734237, t = 0.000048445123372848407605,
x = 0.99980626048671753021, y = 0.99999997568980723193], 1493]


Ссылка к DirectSearch http://forum.exponenta.ru/viewtopic.php ... rectsearch

teolog
Сообщения: 111
Зарегистрирован: Вс янв 22, 2012 8:41 pm

Сообщение teolog » Чт янв 26, 2012 10:52 pm

:o Вы похоже великолепно знаете Мепл. Огромное спасибо! Буду сейчас разбираться с Вашим примером.
И, если не затруднит, могли бы Вы ответить ещё на один вопрос.
Вы мне давали ссылку http://forum.exponenta.ru/viewtopic.php?p=47873&highlight=#47873.
sys:= diff(x1(t),t)=-4*x1(t)^2, diff(x2(t),t)=x1(t)^2-x2(t);
inicond:=x1(0)=1, x2(0)=1;
fens:={x1(t),x2(t)};
sol:=dsolve({sys,inicond},fens,type=numeric);
L:=convert([seq(sol(i),i=0..2.65, 0.05)], matrix);
fd := fopen(`d:\\testfile4.txt`,WRITE);
for i from 1 by 1 to 12 do
for j from 1 by 1 to 3 do
fprintf(fd,"%a %s",L[i,j]," ");
end do;
fprintf(fd,"%s\n","")
end do:
fclose(fd);

Производится запись координат точек в файл (t=, x(t)=, y(t)=).
Я аппроксимирую полученные значения некоторой кривой.
pnts:=[массив значений x],[ массив значений у];
g:=interp(pnts, x);
with(stats):
h:=fit[leastsquare[[x,y], y=a3*x^3+a2*x^2+a1*x+a0]]([pnts]);

В файл записываются три столбца. В переменной pnts - две строчки. Как соединить решение системы дифуров и аппроксимацию? Обязательно ли делать это через файл? можно ли взять координаты только x и y (без времени) и эти два массива поместить в переменную pnts ?
Если Вас не затруднит, могли бы для моего случая показать пример такого склеивания двух кусков кода?

Markiyan Hirnyk
Сообщения: 1366
Зарегистрирован: Вс дек 04, 2011 11:07 pm

Сообщение Markiyan Hirnyk » Чт янв 26, 2012 11:13 pm

teolog писал(а)::o Вы похоже великолепно знаете Мепл. Огромное спасибо! Буду сейчас разбираться с Вашим примером.
И, если не затруднит, могли бы Вы ответить ещё на один вопрос.
Если Вас не затруднит, могли бы для моего случая показать пример такого склеивания двух кусков кода?
Не понял вопрос.

teolog
Сообщения: 111
Зарегистрирован: Вс янв 22, 2012 8:41 pm

Сообщение teolog » Чт янв 26, 2012 11:29 pm

Решая систему ДУ, мы получаем 3 столбца значений: времени, координаты х, и координаты у. И записываем эти три столбца в файл.
Можно ли столбцы х и у поместить соответственно в две строчки автоматически, в пределах одной программы?
pnts:=[массив значений x],[ массив значений у];
В файле, куда помещаются координаты точки на каждом шаге будет такая структура:
t=0, x(t)=1, y(t)=1
t=0.05, x1(t) = .833333228373074 x2(t) = .991809947558016
t = .10 x1(t) = .714285719766109 x2(t) = .972431637562123
t=... , x(t)=... , y(t)=...
И так далее (число строчек равно числу шагов).
А мне надо в итоге получить запись такого вида:
pnts:=[1, .833333228373074, .714285719766109, оставшиеся значения х],[1, .991809947558016, .972431637562123, оставшиеся значения у];
Как это сделать? и можно ли обойтись без записи значений в файл, а сразу же из решения ДУ передать координаты точек в две строчки?

Markiyan Hirnyk
Сообщения: 1366
Зарегистрирован: Вс дек 04, 2011 11:07 pm

Да

Сообщение Markiyan Hirnyk » Пт янв 27, 2012 11:00 am

Это сделать можно:
sys := diff(x1(t), t) = -4*x1(t)^2, diff(x2(t), t) = x1(t)^2-x2(t);
F := dsolve({sys, x1(0) = 1, x2(0) = 1}, {x1(t), x2(t)}, numeric, output = array([0, .25, .5, .75, 1]));
op([1], F);
op([1, 3], F);
op([1, 3, 2], F);
A := LinearAlgebra:-Transpose(convert(rhs(op([1, 3, 2], F)), Matrix));
pnts := convert(A[2, () .. ()], list), convert(A[3, () .. ()], list);

[1., .500000034736135746, .333333310619616519, .249999989309707848, .200000026359693023], [1., .886286436916848386, .726491436776205978, .584006038214527012, .465782044145892416]
Полагаю, что для решения Ваших задач аппроксимация не обязательна.

teolog
Сообщения: 111
Зарегистрирован: Вс янв 22, 2012 8:41 pm

Сообщение teolog » Пт янв 27, 2012 1:10 pm

С вариантом, предложенным Вами, это действительно необязательно. Ваш вариант более точный. Огромное спасибо, Вы мне очень сильно помогли!