MATLAB: ode45 уменьшить шаг по времени

electric_retard · 31/12/13 148

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

GAA · 12/07/07 4529

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

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

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

electric_retard · 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 не исключает вышеописаной проблемы: в данный момент в коде ловится верхняя точка траектории и момент касания земли.
Единственное решение — ограничить шаг по времени на всем диапазоне расчёта, но это увеличивает время расчёта. В окрестности какой-то точки получается невозможно его измельчить?

GAA · 12/07/07 4529

Что-то мешает убрать генерацию ошибки?
Если ничего, то уберите генерацию ошибки (т.е. уберите 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. Но в более сложном варианте их анализ, возможно, понадобится.

electric_retard · 31/12/13 148

GAA
Идею понял, спасибо за подробные разъяснения, так бы дольше на первое сообщение медитировал. Пробую прикрутить.

GAA в сообщении #1101776 писал(а):

Доопределим нулём правую часть уравнения при $x < 0$ .

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

-- 25.02.2016, 16:15 --

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

код: [ скачать ] [ спрятать ]

Используется синтаксис Matlab M

0000000000

0000000000

0000000000                    0.0000000000                    0.0000000000

0000000000

0000401902

0000602853

0001607607

0001786230

0002009509

0002009509

0000003305                    0.0000003305                    0.0002009509

0004019018

0005023773

0010047545

0010940660

0012057054

0012057054

0000016523                    0.0000019827                    0.0010047545

0022104598

0027128369

0052247222

0056712795

0062294760

0062294760

0000082613                    0.0000102440                    0.0050237706

0112532452

0137651264

0263245287

0285573093

0313482843

0313482843

0000413063                    0.0000515503                    0.0251188083

0564670551

0690263562

1318227680

1429865255

1569412016

1569412016

0002065316                    0.0002580819                    0.1255929173

2825331822

3453270649

6592941379

7151092815

7848776896

7848776884

0010326581                    0.0012907400                    0.6279364868

4127907600

7266946137

2961554568

5751401876

9238579529

9238578044

0051632905                    0.0064540305                    3.1389801160

0622527934

6301340106

4680913511

8607542621

6012403545

6012217905

0258164525                    0.0322704830                    15.6773639861

2639891437

0625625881

0206864811

9292188069

5545949924

5522751215

1290822623                    0.1613527453                    77.9510533310

1423083732

1289720223

3717078532

7580649269

0209721727

7316431629

6454113115                    0.8067640568                    379.1793680413

3446280930

8225832057

6555148882

8802476503

5793323451

4535728529

2270565576                    4.0338206144                    1664.7219296900

9538199422

0756903989

3697267978

0971723915

6189409162

0802582112

1200673757                    11.1538879901                   2662.6266853583

9652711803

6559107268

0213636691

5531320163

2190028135

6869136330

9609793593

9112810373

5703071120

8819885856

8816961647

6134058340

2608543401                    19.4147423302                   1898.5331476228

7597032552

3121895578

2383797687

2836422551

9281311347

0878984657

2608543401                    27.6755966703                   979.4744926317

1445930339

4256427592

6963416713

3234266962

1698679553

8031721915

1255376422                   39.8011343125                   146.7152737258

-4.8257370574                   34.9753972551                   116.8605562235

-0.0284745594                   34.9469226957                   0.0041439110

-0.0008657073                   34.9460569884                   0.0000003031

0004261226                    34.9464831110                   0.0000007777

-0.0000000001                   34.9464831109                   0.0000000000

-0.0002130613                   34.9462700497                   -0.0000001642

0001065306                    34.9463765803                   0.0000001383

0001065306                    34.9464831109                   0.0000000260

0000000000                    34.9464831109                   0.0000000000

0000000000                    34.9464831109                   0.0000000000

1111477324

0688789013

3829575771

2792932885

0576715831

2155184975

9475191301                   54.8940022410                   -2035.4523549093

1495763571

5504879568

4563320469

2055390496

9047716445

2599919859

7570499716                   68.6510522126                   -4525.9555265115

-2.6032333465

1862824863

Error using get_env (line 28)

Текущая высота -38073 м выходит за границы диапазона -2:15 км

GAA · 12/07/07 4529

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

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

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

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

electric_retard · 31/12/13 148

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

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

MATLAB: ode45 уменьшить шаг по времени

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