Презентация на тему «Лекция система дифференциальные уравнения»

Физика

Вычислительная

Ю.Н. Прошин кафедра теоретической физики Казанского федерального университета yurii.proshin@kpfu.ru 2004-2013, Казань

Matlab: решение дифференциальных уравнений и другие демонстрации (

Простой пример

Возможный формат вызова процедуры решателя в MatLab: [T,Y]=ode45(@DiffEquatFunc,[Tstart,Tfinal],StartVector).

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

С начальным условием

И аналитическим решением

Простой пример

Снимок экрана, который соответствует численному решению этой задачи в системе MatLab.

Файл-функция, описывающая правую часть уравнения, – текстовый файл с

расширением func1.m – содержит всего две строки function dd=func1(x,y) % название dd=-2*x*y; % правая часть ДУ Знаком % начинаются комментарии. Вызываться такая функция может из другой программы, функции, или, как в этом случае, из командного окна >> [T,Y]=ode45(@func1,[0,2],1); Здесь задан временной интервал от Tstart=0 до Tfinal=2 и начальное значение функции StartVector=1. График полученной таким образом функции Y(T) воспроизводится вызовом встроенной функции plot >> plot(T,Y); Следующей строкой мы кружочками нарисовали на том же графике точное решение в точках полученного вектора-столбца T: >> hold on; plot(T,exp(-T.^2),'ro'); hold off

Заметим, что уравнение (1) можно решить в MatLab и символьно

Приведем часть командного окна, где была вызвана стандартная процедура dsolve >> dsolve('Dy=-2*t*y','y(0)=1') ans = 2 exp(-t ) Здесь также использовано начальное условие.

Видим, что с точностью до переобозначения x ? t результат совпадает с приведенным выше.

Решатели диф

уравнений в MatLab (solvers)

Для решения систем ОДУ в MatLAB реализованы различные методы. Их реализации названы решателями ОДУ. Решатели реализуют следующие методы решения систем дифференциальных уравнений:

Все решатели (ode45, ode23, ode133, ode15s, ode23s, ode23t, ode23tb ) могут решать системы уравнений явного вида y’ = F(t, y). Решатели ode15s, ode23s, ode23t, ode23tb могут решать уравнения неявного вида F(t, y, y’ ) = 0.

Решатели диф

уравнений в MatLab (solvers)

ode45 – одношаговые явные методы Рунге-Кутта 4-го и 5-го порядка {начальная проба решения}. Во многих случаях он дает хорошие результаты; ode23 – одношаговые явные методы Рунге-Кутта 2-го и 3-го порядка. При умеренной жесткости системы ОДУ и низких требованиях к точности этот метод может дать выигрыш в скорости решения; ode133 – многошаговый метод Адамса-Башворта-Мултона переменного порядка. Это адаптивный метод, который может обеспечить высокую точность решения; ode15s – многошаговый метод переменного порядка (от 1-го до 5-го, по умолчанию 5), использующий формулы численного дифференцирования. Это адаптивный метод, его стоит применять, если решатель ode45 не обеспечивает решения; ode23s – одношаговый метод, использующий модифицированную формулу Розенброка 2-го порядка. Может обеспечить высокую скорость вычислений при низкой точности; ode23t – метод трапеций с интерполяцией. Этот метод дает хорошие результаты при решении задач, описывающих осцилляторы с почти гармоническим выходным сигналом; ode23tb – неявный метод Рунге-Кутта в начале решения и метод, использующий формулы обратного дифференцирования 2-го порядка в последующем. При низкой точности этот метод может оказаться более эффективным, чем ode15s.

Движение заряженной частицы

Закон Кулона

Пусть некоторая точка массы m с зарядом q движется в электрическом поле двух неподвижных зарядов Q1 и Q2

Движение заряженной частицы

Закон Кулона

Движение заряженной частицы

Закон Кулона

Пусть масса частицы m = 1, ее заряд q = 1. Перейдем к безразмерным единицам, и будем считать, что данная задача является "плоской". Введем следующие обозначения:

,

,

И

.

Движение заряженной частицы

Закон Кулона

Пусть масса частицы m = 1, ее заряд q = 1. Перейдем к безразмерным единицам, и будем считать, что данная задача является "плоской". Введем следующие обозначения:

,

,

И

.

Движение заряженной частицы

Закон Кулона

Рассмотрим простейший случай финитного движения с Q1 = –50, Q2 = 0, С1 = (5,0) и С2 = (0,10). При таких начальных параметрах (Q2 = 0 и Q1 < 0) наша точка движется в притягивающем поле только первого заряда и, как мы помним из классической механики, должна описывать вокруг него эллипс. Проверим, запишем правую часть системы уравнений как файл-функцию, назвав ее pointq12.

function f=pointq12(t,x) global Q1 Q2 C1x C1y C2x C2y f=[x(3);x(4);... Q1*(x(1)-C1x)/(sqrt((x(1)-C1x)^2+(x(2)-C1y)^2))^3+... Q2*(x(1)-C2x)/(sqrt((x(1)-C2x)^2+(x(2)-C2y)^2))^3;... +Q1*(x(2)-C1y)/(sqrt((x(1)-C1x)^2+(x(2)-C1y)^2))^3+... Q2*(x(2)-C2y)/(sqrt((x(1)-C2x)^2+(x(2)-C2y)^2))^3];

И

.

Движение заряженной частицы

Закон Кулона

Решим систему дифференциальных уравнений, вызвав процедуру ode45 из "пустой" файла-функции pointDyn.m

Function pointdyn() clear all global Q1 Q2 c1x c1y c2x c2y Q1=-50; Q2=-0.; C1x=5; c1y=0; c2x=0; c2y=10; x0=0; y0=0; vx0=0; vy0=4.3; T1=4000; [t,h]=ode45(@pointq12,[0,t1],[x0,y0,vx0,vy0]); x=h(:,1); y=h(:,2); x1=c1x; y1=c1y; x2=c2x; y2=c2y; plot(x,y,'b-'); % отрисовка траектории hold on % отрисовка положения неподвижных зарядов plot(x1,y1,'r+',x2,y2,'r*','markersize',15); plot(x1,y1,'ro',x2,y2,'ro','markersize',15); % comet(x,y); % отрисовка "движения" hold off;

И

.

Движение заряженной частицы

Закон Кулона

Для данного примера относительной точности 10-3, заложенной по умолчанию в процедуре ode45, недостаточно. Придется либо уменьшать этот параметр, либо пробовать другие процедуры.

И

.

Движение заряженной частицы

Закон Кулона

Изменим относительную точность решения на три порядка ? 10-6 : tol = 1e-6; [t,h]=ode45(@pointq12,[0,T1],[x0,y0,vx0,vy0],...odeset('RelTol',tol))

И

.

Движение заряженной частицы

Закон Кулона

Время расчета увеличилось, но теперь мы получили вполне приемлемый результат. Теперь можно поэкспериментировать с начальными условиями и зарядами тел (например, можно убедиться, что при последовательном увеличении на единицу заряда Q1 с –50 до –46 движение становится инфинитным). Естественно, что движение станет также инфинитным, если взять заряды одного знака, т.е. Q1 > 0.

И

.

Движение заряженной частицы

Закон Кулона

Введем значение Q2 = –0.2. Это внесет возмущение в орбиту движущейся точки.

И

.

Движение заряженной частицы

Закон Кулона

Введем значение Q2 = +0.2. Траектория становится незамкнутой .

И

.

Движение заряженной частицы

Закон Кулона

Q1 = –50 и Q2 = –1.5. T1 = 8000.

.

Движение заряженной частицы

Закон Кулона

Следует подчеркнуть, что мы использовали модель точечных зарядов, т.е. пренебрегали возможностью "попадания" зарядов друг в друга. Дальнейшее улучшение программы связано с контролем в ходе решения выполнения закона сохранения энергии (особенно это существенно при решении задачи многих тел). Теперь можно экспериментировать самостоятельно ! Ю.Н. Прошин, И.М. Еремин. Вычислительная физика (Практический курс) Казань: Казанский государственный университет, 2009. ("Задание 4. Движение заряда в кулоновском поле")

И

.

Движение под действием сил тяжести и трения

Рассмотрим траекторию движения пули под действием силы тяжести. При отсутствии сопротивления воздуха это будет парабола. При скорости пули больше скорости звука сила сопротивления воздуха пропорциональна квадрату скорости и противоположна направлению движения. Уравнение движения пули массой m будет следующим

Примем для простоты, что коэффициент пропорциональности k в силе трения зависит от плотности воздуха ?, которая, в общем случае, может меняться с высотой y, площади поперечного сечения пули S и некоторого постоянного безразмерного параметра b порядка единицы, учитывающего форму пули. Из соображений размерности k = b?S.

И

.

Движение под действием сил тяжести и трения

Пусть масса пули - m = 9 грамм, S = 0.5 см2 (~ калибр 7.62 мм). Пусть ? = ?(0) = 1.22 кг/м3, g = 9.8 м/с2, коэффициент b = 0.5. При t0 = 0: x0=0, y0=0, а vx0 = 800 м/с, vy0 = 100 м/с. Переведем все в систему Си и обезразмерим.

По осям - метры

- - - - ode45, + + + + ode113.

И

.

Движение под действием сил тяжести и трения

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

И

.

Движение под действием сил тяжести и трения

I. Как изменится траектория пули с учетом распределения плотности воздуха по высоте. Построить аппроксимацию ?(y) по следующим данным: в Европе плотность воздуха у поверхности Земли равна 1.258 кг/м3, на высоте 5 км – 0.735 кг/м3, на высоте 20 км – 0.087 кг/м3, на высоте 40 км – 0.004 кг/м3 II. При каком угле вылета пуля достигает максимальной дальности? III. Если одну пулю выстрелить горизонтально из ствола, а другую пулю бросить с той же самой высоты в тот же самый момент, упадут ли обе из них в одно и то же время? IV. Если пулю выстрелить из винтовки вертикально вверх, какой будет ее окончательная скорость, когда она попадет в макушку чьей-то головы во время своего полета вниз? Построить фазовую диаграмму (y, Vy) при разных параметрах задачи. Зависит ли дальность от массы пули? Какая пуля улетит дальше, более тяжелая или более легкая? Учесть осевое вращение пули, возникающее в нарезном стрелковом оружии. Как это влияет на дальность и точность стрельбы?

Теперь можно проводить компьютерные "эксперименты", меняя параметры задачи. Например, можно учесть изменение плотности воздуха с высотой, или даже осевое вращение пули, возникающее в нарезном стрелковом оружии и оценить как это влияет на точность стрельбы (см. Вычислительная физика "Задание 2. Полет пули", на стр. 26).

И

.

Осциллятор Ван дер Поля

Жесткость системы ОДУ

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

здесь ? – параметр. Перепишем

И

.

Осциллятор Ван дер Поля

Жесткость системы ОДУ

Используем в качестве функции, описывающей систему ДУ (13), функцию vanderpoldemo, входящую в стандартный демонстрационный пример MatLab – odedemo.

function dydt = vanderpoldemo(t,y,Mu) %VANDERPOLDEMO Defines the van der Pol equation % Copyright 1984-2002 The MathWorks, Inc. % $Revision: 1.2 $ $Date: 2002/06/17 13:20:38 $ dydt = [y(2); Mu*(1-y(1)^2)*y(2)-y(1)];

И

.

Осциллятор Ван дер Поля

Жесткость системы ОДУ

Для малых ? порядка единицы практически любой MatLab решатель ОДУ сможет эффективно решить уравнение Ван дер Поля. Для больших значений ? > 100 система ОДУ становится жесткой. Для быстрого и эффективного интегрирования таких систем должны быть использованы специальные методы, реализованные в ode15s, ode23s, ode23t, ode23tb. Сравним работу двух процедур ode45 (синяя сплошная линия на графиках) и ode15s (прерывистая красная) при разных значениях ?… Начальные условия

И

.

Осциллятор Ван дер Поля

Жесткость системы ОДУ

tspan = [0, 100]; x0 = [0.5; 0]; Mu = 0.0; disp(['Fig0 tspan = [0, 100]; mu=', num2str(Mu)]) tic % Засекаем время [t,x] = ode45(@vanderpoldemo, tspan, x0,[],Mu); toc % Останавливаем и печатаем время % Plot of the solution plot(t,x(:,1),'b','LineWidth',4) xlabel('t'); ylabel('solution x') title(['van der Pol Equation, \mu = ', num2str(Mu)]) hold on; tic [t,x] = ode15s(@vanderpoldemo, tspan, x0,[],Mu); toc; plot(t,x(:,1),'r--','LineWidth',4); hold off

И

.

Осциллятор Ван дер Поля

Жесткость системы ОДУ

tspan = [0,100]; mu=0 ode45 => 0.103 sec. ode15s => 0.284 sec.

Периодическое решение для гармонического осциллятора с амплитудой 0.5

И

.

Осциллятор Ван дер Поля

Жесткость системы ОДУ

tspan = [0,100]; mu=0.1 ode45 => 0.112 sec. ode15s => 0.447 sec.

И

.

Осциллятор Ван дер Поля

Жесткость системы ОДУ

tspan=[0,100]; mu=1 ode45 => 0.191 sec. ode15s => 0.796 sec.

И

.

Осциллятор Ван дер Поля

Жесткость системы ОДУ

Система становится жестче – малое изменение параметра, приводит к сильному изменению функции. Хотя по-прежнему ode45 быстрее ode15s.

tspan=[0,100]; mu=10 ode45 =>0.450 sec. ode15s => 1.07 sec.

И

.

Осциллятор Ван дер Поля

Жесткость системы ОДУ

tspan =[0,2000]; mu = 175 ode45 => 190.4 sec. ode15s => 2.631 sec.

И

.

Осциллятор Ван дер Поля

Жесткость системы ОДУ

При ? = 1000 ode45 отказывается работать, а время работы ode15s увеличился на треть при увеличении временного диапазона почти в 5 раз:

tspan=[0,10000]; mu=1000 ode45 => Not solved! ode15s > 3.166 sec.

И

.

Осциллятор Ван дер Поля

Жесткость системы ОДУ

Причина в том, что при увеличении параметра ? начинают сильно различаться порядки коэффициентов при разных слагаемых. Именно степень этого различия чаще всего и определяет жесткость системы ОДУ

В качестве соответствующей характеристики выбирают матрицу Якоби (якобиан) векторной функции F(t,X), определяющей правую часть системы ОДУ. Чем сильнее вырождена матрица Якоби, т.е. функциональная матрица, составленная из производных F(t,X), тем жестче система уравнений.

При ? = 1000 ode45 отказывается работать, а время работы ode15s увеличился на треть при увеличении временного диапазона почти в 5 раз:

tspan=[0,10000]; mu=1000 ode45 => Not solved! ode15s > 3.166 sec.

И

.

Система Лоренца (С

П.Кузнецов, Динамический хаос)

Система Лоренца

Система Лоренца (Г

Шустер, Детерминированный хаос)

Система Лоренца (Г

Шустер, Детерминированный хаос)

Система Лоренца (Г

Шустер, Детерминированный хаос)

Система Лоренца (Г

Шустер, Детерминированный хаос)

Система Лоренца (Г

Шустер, Детерминированный хаос)

Система Лоренца (Г

Шустер, Детерминированный хаос)

Система Лоренца (Г

Шустер, Детерминированный хаос)

Система Лоренца (Г

Шустер, Детерминированный хаос)

Система Лоренца (Г

Шустер, Детерминированный хаос)

Система Лоренца (Г

Шустер, Детерминированный хаос)

Система Лоренца (Г

Шустер, Детерминированный хаос)

Система Лоренца (Г

Шустер, Детерминированный хаос)

Система Лоренца (C

Кузнецов, Динамический хаос)

Система Лоренца (C

Кузнецов, Динамический хаос)

Система Лоренца (C

Кузнецов, Динамический хаос)

Система Лоренца (C

Кузнецов, Динамический хаос)

Система Лоренца (C

Кузнецов, Динамический хаос)

Система Лоренца (C

Кузнецов, Динамический хаос)

Система Лоренца (C

Кузнецов, Динамический хаос)

Система Лоренца (C

Кузнецов, Динамический хаос)

O

O1

O2

Система Лоренца

function lorenz(action) %LORENZ Plot the orbit around the Lorenz chaotic attractor. % This demo animates the integration of the three coupled nonlinear differential % equations that define the "Lorenz Attractor", a chaotic system first described % by Edward Lorenz of the Massachusetts Institute of Technology. % Copyright 1984-2005 The MathWorks, Inc. % $Revision: 5.13.4.3 $ $Date: 2005/12/15 20:52:53 $ % The values of the global parameters are global SIGMA RHO BETA SIGMA = 10.; BETA = 8./3.; RHO = 30; % % 28 1 13.927 24.06 24.74 …………………………………………………

Система Лоренца

? = SIGMA = 10; b = BETA = 8/3 = 2.66; r = RHO = 30

function lorenz(action) %LORENZ Plot the orbit around the Lorenz chaotic attractor. % This demo animates the integration of the three coupled nonlinear differential % equations that define the "Lorenz Attractor", a chaotic system first described % by Edward Lorenz of the Massachusetts Institute of Technology. % Copyright 1984-2005 The MathWorks, Inc. % $Revision: 5.13.4.3 $ $Date: 2005/12/15 20:52:53 $ % The values of the global parameters are global SIGMA RHO BETA SIGMA = 10.; BETA = 8./3.; RHO = 30; % % 28 1 13.927 24.06 24.74 …………………………………………………

Система Лоренца

RHO = 0.9 < 1

O

Система Лоренца

RHO = 7.5 < 13.927

O1

Система Лоренца

RHO = 7.5 < 13.927

O2

Система Лоренца

RHO = 19 > 13.927

O2

Система Лоренца

24.06 < RHO = 24.5 < 24.74

O2

Система Лоренца

24.06 < RHO = 24.5 < 24.74

O1

Система Лоренца

24.06 < RHO = 24.5 < 24.74

? Множество ?1 ? странный аттрактор

Система Лоренца

RHO = 28 > 24.74

Система Лоренца

Система Лоренца

Система Лоренца

Литература

Ю.Н. Прошин, И.М. Еремин. Вычислительная физика (Практический курс) - Казань: Казанский государственный университет, 2009. – 180 с. С. П. Кузнецов, Динамический хаос.- М: ФИЗМАТЛИТ, 2006. - 356 с. А.Ю. Лоскутов, А.С. Михайлов, Основы теории сложных систем. – М.–Ижевск: Институт компьютерных исследований, 2007. – 620 с. Г. Шустер Детерминированный Хаос: Введение. -М.: Мир, 1988. - 253 с.

The End