2014 dxdy logo

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

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




Начать новую тему Ответить на тему
 
 MATLAB: ode45 уменьшить шаг по времени
Сообщение19.02.2016, 16:44 
Аватара пользователя


31/12/13
148
Решаются уравнения внешней баллистики, в процессе решения используется кусочно-гладкая функция, заданная не некотором ограниченном интервале. ode45 так и норовит вылезти за границы этого диапазона. При решении ОДУ в MATLAB есть возможность ограничения максимального шага по времени, но мне хотелось бы иметь возможность ограничивать шаг по времени только, когда одна из переменных приближается к некоторому значению, а не всё время решения ДУ. Я хочу невозможного? :-)

 Профиль  
                  
 
 Re: MATLAB: ode45 уменьшить шаг по времени
Сообщение19.02.2016, 20:33 
Заслуженный участник


12/07/07
4438
Проблема, на мой взгляд, описана не достаточно ясно — крайне желательно формулировать подробней. Но на первый взгляд все выглядит просто.

При «пересечении» указанной переменной заданного значения завершать интегрирование с первоначально выбранным шагом по времени и стартовать продолжение расчета с меньшим шагом по времени. Для отслеживания переключения на другой шаг по времени использовать events.

Чем не годится такой стандартный вариант?

 Профиль  
                  
 
 Re: MATLAB: ode45 уменьшить шаг по времени
Сообщение24.02.2016, 11:23 
Аватара пользователя


31/12/13
148
GAA
На каждом шаге используется зависимость температуры от высоты:
Используется синтаксис Matlab M
if H >= -2000 && H <= 9324
    tau_N      = 289 - 6.328e-3*h;
   
elseif H > 9324 && H <= 12e3
    tau_N      = 230 - 6.328e-3*(H-9324)+1.172e-6*(H-9324).^2;
   
elseif H > 12e3 && H <= 15e3
    tau_N      = 221.5;
else
    error('Текущая высота %5.0f м выходит за границы диапазона -2:15 км',H);
end

Соответственно, если отпускаем шаг по времени в свободное плавание, то можно достаточно часто ловить ошибку, т.к. сама задача допускает большой шаг. Пример: навесная траектория, скорость при касании 600 м/с, максимальный шаг по времени 12 с. За это время тело пролетит много больше 2000 м. Сам решатель допускает еще больший шаг по времени, 12 с — это ограничение, введеное мной (1/10 часть от времени движения тела по параболе по школьной формуле).
Создание events не исключает вышеописаной проблемы: в данный момент в коде ловится верхняя точка траектории и момент касания земли.
Единственное решение — ограничить шаг по времени на всем диапазоне расчёта, но это увеличивает время расчёта. В окрестности какой-то точки получается невозможно его измельчить?

 Профиль  
                  
 
 Re: MATLAB: ode45 уменьшить шаг по времени
Сообщение24.02.2016, 16:29 
Заслуженный участник


12/07/07
4438
Что-то мешает убрать генерацию ошибки?
Если ничего, то уберите генерацию ошибки (т.е. уберите error(...)), если она стоит в одной из вызываемых ode45 функций. Вместо этого доопределите функцию за границами.

1. Если нужно менять максимальный шаг по времени в зависимости от значения момента времени, то это делается очевидным разбиением всего промежутка интегрирования [0, Tf] на отрезки с заданием на каждом отрезке своего максимального шага по времени.

2. Если разные максимальные шаги нужно выбирать в зависимости от значения решения, то можно использовать events.

Не очень уместный пример. (Более уместный с ходу не придумал).
Пусть надо решить задачу $x’ = -\sqrt{x}$, $x(0) =1$. При этом хотим, чтобы шаги по времени убывали по мере приближения $x$ к нулю.
Доопределим нулём правую часть уравнения при $x < 0$.
код: [ скачать ] [ спрятать ]
Используется синтаксис Matlab M
function EventsDemo
  global h
  tau = 1; h = 0.5;
  opts=odeset('MaxStep', tau, 'Events', @eventsfun);    
  T = 0; Tf = 2.5; X = 1;
  Ts = [T];
  Xs = [X];
  while X(end) > 1E-16
    [T, X, TE, XE, IE] = ode45(@f, [T(end), Tf], X(end), opts);
    if length(T) > 1
       Ts = [Ts; T(2:end)];
       Xs = [Xs; X(2:end)];
       tau = tau /2; h = h - 0.1;        
       opts=odeset('MaxStep', tau, 'Events', @eventsfun);
    end
  end
  plot(Ts, Xs, 'ok');
end

function y= f(t, x)
  if x > 0
     y = -sqrt(x);
  else
     y = 0;
  end  
end
   
function [value, isterminal, direction] = eventsfun(t, x)
 global h
  value = x-h;
  isterminal = 1;    
  direction = 0;
end
Я не анализировал значения TE, XE, IE. Но в более сложном варианте их анализ, возможно, понадобится.


Вложения:
Комментарий к файлу: Полученный в Matlab 6.5 график решения при выполнении приведенной в сообщении функции.
plot.png [10.58 Кб]
Скачиваний: 0
 Профиль  
                  
 
 Re: MATLAB: ode45 уменьшить шаг по времени
Сообщение25.02.2016, 15:27 
Аватара пользователя


31/12/13
148
GAA
Идею понял, спасибо за подробные разъяснения, так бы дольше на первое сообщение медитировал. Пробую прикрутить.
GAA в сообщении #1101776 писал(а):
Доопределим нулём правую часть уравнения при $x < 0$.

Вот видите, в Вашем случае вообще как шаг не мельчи, без доопределения ерунда получится (ну или почти ерунда), ведь решатель будет смотреть значения функции и за двойкой.
Функция в моем листинге выдаёт температуру атмосферы в зависимости от высоты над пов-тью земли, взята из книги. Для высот выше можно доопределить, а вот температура атмосферы даже на -2000 м — это что-то такое не очень понятное, слабо представляю как ее там доопределить и в то же время от нее будет зависеть результат.

-- 25.02.2016, 16:15 --

Вот кстати типичный вывод по ходу решения:
Строки, где столбец один — это вызовы функции параметров атмосферы, показана текущая высота.
Где три столбца — это вызов event функции, первый столбец — разница во врмени между вызовами event функции, второй — текущее время, третий — разница в высоте между вызовами.
код: [ скачать ] [ спрятать ]
Используется синтаксис Matlab M
0.0000000000
0.0000000000

0.0000000000                    0.0000000000                    0.0000000000

0.0000000000
0.0000401902
0.0000602853
0.0001607607
0.0001786230
0.0002009509
0.0002009509

0.0000003305                    0.0000003305                    0.0002009509

0.0004019018
0.0005023773
0.0010047545
0.0010940660
0.0012057054
0.0012057054

0.0000016523                    0.0000019827                    0.0010047545

0.0022104598
0.0027128369
0.0052247222
0.0056712795
0.0062294760
0.0062294760

0.0000082613                    0.0000102440                    0.0050237706

0.0112532452
0.0137651264
0.0263245287
0.0285573093
0.0313482843
0.0313482843

0.0000413063                    0.0000515503                    0.0251188083

0.0564670551
0.0690263562
0.1318227680
0.1429865255
0.1569412016
0.1569412016

0.0002065316                    0.0002580819                    0.1255929173

0.2825331822
0.3453270649
0.6592941379
0.7151092815
0.7848776896
0.7848776884

0.0010326581                    0.0012907400                    0.6279364868

1.4127907600
1.7266946137
3.2961554568
3.5751401876
3.9238579529
3.9238578044

0.0051632905                    0.0064540305                    3.1389801160

7.0622527934
8.6301340106
16.4680913511
17.8607542621
19.6012403545
19.6012217905

0.0258164525                    0.0322704830                    15.6773639861

35.2639891437
43.0625625881
82.0206864811
88.9292188069
97.5545949924
97.5522751215

0.1290822623                    0.1613527453                    77.9510533310

175.1423083732
213.1289720223
402.3717078532
435.7580649269
477.0209721727
476.7316431629

0.6454113115                    0.8067640568                    379.1793680413

847.3446280930
1013.8225832057
1847.6555148882
2009.8802476503
2176.5793323451
2141.4535728529

3.2270565576                    4.0338206144                    1664.7219296900

2800.9538199422
3063.0756903989
4433.3697267978
4799.0971723915
5038.6189409162
4804.0802582112

7.1200673757                    11.1538879901                   2662.6266853583

5434.9652711803
5656.6559107268
6832.0213636691
7193.5531320163
7364.2190028135
7061.6869136330
5295.9609793593
5484.9112810373
6448.5703071120
6682.8819885856
6840.8816961647
6702.6134058340

8.2608543401                    19.4147423302                   1898.5331476228

6981.7597032552
7080.3121895578
7566.2383797687
7672.2836422551
7742.9281311347
7682.0878984657

8.2608543401                    27.6755966703                   979.4744926317

7859.1445930339
7879.4256427592
7932.6963416713
7937.3234266962
7880.1698679553
7828.8031721915

12.1255376422                   39.8011343125                   146.7152737258


-4.8257370574                   34.9753972551                   116.8605562235


-0.0284745594                   34.9469226957                   0.0041439110


-0.0008657073                   34.9460569884                   0.0000003031


0.0004261226                    34.9464831110                   0.0000007777


-0.0000000001                   34.9464831109                   0.0000000000


-0.0002130613                   34.9462700497                   -0.0000001642


0.0001065306                    34.9463765803                   0.0000001383


0.0001065306                    34.9464831109                   0.0000000260


0.0000000000                    34.9464831109                   0.0000000000


0.0000000000                    34.9464831109                   0.0000000000

7683.1111477324
7507.0688789013
6480.3829575771
6209.2792932885
5872.0576715831
5910.2155184975

19.9475191301                   54.8940022410                   -2035.4523549093

5326.1495763571
4920.5504879568
2478.4563320469
1707.2055390496
1029.9047716445
1384.2599919859

13.7570499716                   68.6510522126                   -4525.9555265115

-2.6032333465
864.1862824863
Error using get_env (line 28)
Текущая высота -38073 м выходит за границы диапазона -2:15 км

 Профиль  
                  
 
 Re: MATLAB: ode45 уменьшить шаг по времени
Сообщение28.02.2016, 21:48 
Заслуженный участник


12/07/07
4438
Если решение не выходит за границы промежутка, на котором задана Ваша функция, то всё не так страшно, при некритическом рассмотрении. Солвер смотрит за границы только для выбора шага интегрирования, т.е. возвращаемое ode45 решение может зависеть от значений за границами через выбор шага интегрирования.

Проблема в другом. Функция ode45 возвращает решение с заданной точностью, если начальное условие задано точно. Если прерывать решение задачи и стартовать с последней посчитанной точки, то уже ode45 не вернёт решение на всём промежутке интегрирования с заданной точностью.

Т.е. тот вариант, который я предложил изначально, есть не более чем предварительная прикидка решения.

Если же нужно серьёзно решать и тонко рулить шагом интегрирования, то придётся или искать другой (сторонний) солвер, или писать код солвера самостоятельно.

 Профиль  
                  
 
 Re: MATLAB: ode45 уменьшить шаг по времени
Сообщение05.03.2016, 11:19 
Аватара пользователя


31/12/13
148
GAA
Поэкспериментировав я понял, что ode45 на автомате в принципе нельзя разрешать рулить шагом по времени. Изначально считал, что можно, что машина сама адекватно его выставит. В совокупности с заведомо большим Tf и получалась ерунда, т.е. результаты получаются не совсем адекватные при внимательном рассмотрении.
Теперь выбор шага зависит от начальной скорости, в такой постановке никаких аномалий не выскакивает.

 Профиль  
                  
Показать сообщения за:  Поле сортировки  
Начать новую тему Ответить на тему  [ Сообщений: 7 ] 

Модераторы: Karan, Toucan, PAV, maxal, Супермодераторы



Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей


Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете добавлять вложения

Найти:
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group