Федеральное агентство по образованию 3 страница

Математическое моделирование

Математическое моделирование позволяет глубже понять физический процесс, а также рассмотреть условия, которые преимущественно не могут быть достигнуты в эксперименте.

Этапы построения модели:

1. Создаётся физическая модель – осуществляется выбор закона, описывающего явление, оцениваются значения параметров задачи и т. д.

1. Создаётся математическая модель – производится запись уравнений, которые затем требуется решить.
2. Осуществляется выбор алгоритма для решения математической модели.
3. Производится подгонка параметров модели и проверка её адекватности.

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

Пример 1. Радиоактивный распад описывается дифференциальным уравнением

(1)

где - константа распада, N - количество радиоактивных ядер. Решением этого уравнения является выражение

, (2)

из которого видно, что время релаксации - время, за которое количество ядер уменьшается в e раз. Оно связано с периодом полураспада T формулой , которую можно получить из соотношения

Пример 2. Процесс теплопроводности от тела с температурой T в окружающую среду с температурой через материал с теплопроводностью описывается уравнением

. (3)

Как и в предыдущей задаче, разность температур спадает по экспоненте

. (4)

Пример 3. Как известно, при движении в вязкой среде возникает вязкое трение. Уравнение движения записывается в данном случае в виде, аналогичном (1)

, (5)

где – коэффициент, зависящий от вязкости среды, формы тела и его массы.

Соответственно, скорость тела в вязкой среде должна спадать по экспоненте: , где - время релаксации. Заметим, что такая зависимость наблюдается только при ламинарном обтекании тела, при возникновении турбулентных вихрей сила сопротивления становится пропорциональна квадрату скорости

Пример 4. Диффузия также описывается дифференциальным уравнением, аналогичным (3). Рассмотрим процесс диффузии с поверхности

. (6)

Здесь N - количество молекул, продиффундировавших через площадку S за счёт градиента концентрации , - коэффициент диффузии, – характерный масштаб изменения концентрации вещества. В результате диффузии разность концентраций будет спадать по экспоненте.

где - время релаксации.

Пример 5. Разрядка конденсатора с зарядом и ёмкостью через резистор с сопротивлением описывается аналогичным уравнением. Действительно, согласно 2-му правилу Кирхгофа напряжение в такой схеме складывается из напряжений на конденсаторе и на резисторе

Отсюда получим дифференциальное уравнение

. (7)

с решением

.

Пример 6. Спад тока, протекающего через индуктивность с сопротивлением , описывается уравнением того же типа. ЭДС самоиндукции в катушке компенсируется падением напряжения на сопротивлении: , откуда получаем и

. (8)

Решение, как и в предыдущих случаях – экспонента.

Пример 7. Кинетика химических реакций типа описывается теми же уравнениями: изменение концентрации вещества А прямо пропорционально произведению концентраций обоих веществ.

(9)

Отличие состоит в том, что время релаксации здесь определяется концентрацией второго вещества .

Пример 8. Аналогичные модели используются в популяционных задачах. Рассмотрим модель Мальтуса. Численность популяции за время изменяется благодаря рождаемости и смертности . При сложении этих изменений получим дифференциальное уравнение первого порядка

(9)

с решением . В случае, если рождаемость превышает смертность, получаем экспоненциальный рост численности, в обратном случае - спад. При равенстве этих показателей численность не меняется.

Проведём моделирование этой задачи в среде MATLAB. В отдельном файле maltus.m зададим функцию, описывающую дифференциальное уравнение

function dn = maltus (t,n,alpha,beta)

dn = (alpha-beta)*n;

Теперь зададим коэффициенты и используем эту функцию в операторе решения обыкновенных дифференциальных уравнений по схеме Рунге-Кутта 2-3 порядка

alpha = 0.01; beta = 0.02;

[t,n]=ode23(@maltus,[0 20],100,[],alpha,beta);

Рис. 1. Зависимости численности популяции от времени в модели Мальтуса

Результат отобразим в виде графика, подписав оси и дав заголовок графику.

plot(t,n)

xlabel(‘Время ‘)

ylabel(‘Численность’)

title(‘Модель Мальтуса: dN/dt=\gammaN’).

На рисунке 1 показаны результаты расчёта для трёх ситуаций: , , . Как видно из него, за 10 лет численность популяции меняется в e раз.

Нелинейные математические модели

Реальные процессы, как правило, оказываются линейными в очень ограниченной области. Решение нелинейных уравнений аналитическими методами наталкивается на большие трудности, поэтому они решаются, как правило, численными методами.

Задача 1. Популяционная задача с учетом ограничения по ресурсам

В частности, в модели Мальтуса мы получаем неограниченный рост численности популяции, чего в действительности быть не может – численность, как правило, ограничена наличием тех или иных ресурсов, источников питания и т.п. Эта особенность учтена в так называемой логистической модели, описываемой уравнением

(10)

где - предельная численность популяции, которую может прокормить данная территория

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

1. N<Np – идёт экспоненциальный рост, замедляющийся по мере приближения к предельной численности.

2. N=Np – численность не меняется.

3. N>Np - идёт экспоненциальный спад, который замедляется по мере приближения к предельной численности.

Рис. 2. Зависимости численности популяции от времени в логистической модели

Проведём моделирование этой задачи в среде MATLAB. В отдельном файле logistic_model.m зададим функцию, описывающую дифференциальное уравнение

function dn=logistic_model(t,n,gamma,Np)

dn=gamma*(1-n./Np).*n;

Произведём расчёт для трёх различных начальных значений численности, результаты расчёта отобразим на графике (см. рис.2), подпишем оси, вставим заголовок и легенду

[t,y]=ode45(@logistic_model, [0 20], 100,[],0.1, 200);

[t1,y1]=ode45(@logistic_model, [0 20], 200,[],0.1, 200);

[t0,y0]=ode45(@logistic_model, [0 20], 400,[],0.1, 200);

plot(t,y,’:’,t1,y1,’-‘,t0,y0,’-.’),

grid, xlabel('Время t'), ylabel('Численность N')

title('Логистическая модель: dN/dt=\gamma_0\cdot(1-N/N_p)N; \gamma_0=0.1'),

legend(‘N_0<N_p’,’N_0=N_p’,N_0>N_p’)

Как видно из рисунка 2, независимо от начальной численности с течением времени количество особей в популяции стремится к предельному значению

Задача 2. Популяционная задача с учетом ограничения по ресурсам и модуляции параметров

Преимуществом численного моделирования является относительная простота, с которой получаются решения для довольно сложных моделей. Например, можно ввести периодическую вариацию параметров задачи. Рассмотрим случай периодической изменчивости коэффициента прироста . Функция, определяющая дифференциальное уравнение, в этом случае записывается

function dn=logistic_model1(t,n,gamma0,Np,T)

gamma=gamma0*(1+cos(2*pi*t/T);

dn=gamma*(1-n./Np).*n;

Результат расчёта для трёх различных начальных состояний при T=1 показан на рис.3. Как видно из сравнения с рис.2, на кривых появляются модуляции.

Рис. 3. Зависимости численности популяции от времени в логистической модели с модуляцией коэффициента прироста

Рассмотрим случай периодической изменчивости оптимальной численности популяции . Функция, определяющая дифференциальное уравнение, в этом случае записывается

function dn=logistic_model2(t,n,gamma,Np0,T)

Np=Np0*(1+0.2*cos(2*pi*t/T);

dn=gamma*(1-n./Np).*n;

Как видно из рис.4, во всех решениях, независимо от начального состояния, появляется периодичность. Интересно, что несмотря на симметричность модуляции наблюдается снижение оптимального уровня, происходящее, видимо, за счёт нелинейности задачи.

Рис. 4. Зависимости численности популяции от времени в логистической модели с модуляцией оптимального уровня

Задача 3. Нелинейная модель динамики численности популяции

Ещё один фактор, который необходимо учитывать в популяционной динамике, состоит в наличии пола. В таком случае прирост популяции определяется вероятностью встреч особей противоположного пола и поэтому пропорционален квадрату численности популяции (при равенстве долей обоих полов в ней): . Количество смертей, между тем, пропорционально численности популяции: . В результате получим уравнение

(11)

Как видим, уравнение подобно модели Мальтуса с показателем . Его значение зависит от численности популяции, соответственно, возможно три поведения системы:

1. – происходит убыль численности;

2. – состояние не меняется;

3. – идёт быстрый рост популяции.

Алгоритм

Создадим отдельный файл model_nl.m, описывающий функцию для решения дифференциального уравнения

function dn = model _nl(t,n,alpha,beta)

dn = alpha* n.* n – beta * n;

Модели на основе систем обыкновенных дифференциальных уравнений

Задача 1. Популяционная задача с учетом полового состава

Динамика численности популяции с учётом полового состава определяется двумя уравнениями типа (11), образующими систему

(12)

Алгоритм

Создадим отдельный М-файл, где зададим функцию, описывающую ОДУ

function dy=popul(t,y,alpha,beta);

dy=[alpha*y(1)*y(2)-beta*y(1);

alpha*y(1)*y(2)-beta*y(2)];

Получим решение при помощи оператора решателя ОДУ по схеме Рунге-Кутты 4-5 порядков[1] для различных начальных значений численности особей одного пола.

[t,y]=ode45(@popul, [0 20], [100 100],[],0.001,0.1);

[t1,y1]=ode45(@popul, [0 20], [110 100],[],0.001,0.1);

[t0,y0]=ode45(@popul, [0 20], [90 100],[],0.001,0.1);

Графически представим результаты расчетов, разбив графическое окно на два подокна: сверху – зависимость от времени численности особей одного поля, снизу – другого. Подпишем оси, нанесём сетку, в заголовок выпишем значения параметров.

subplot(211), plot(t,y(:,1),t1,y1(:,1),’--‘,t0,y0(:,1),’-.’), grid, xlabel('t'), ylabel('n_1')

title('dn_1/dt=\alphan_1\cdotn_2-\betan_1;\alpha=0.001; \beta=0.1; n_{10}=100; 110; 90; n_{20}=100')

subplot(212), plot(t,y(:,2),t1,y1(:,2),’--‘,t0,y0(:,2),’-.’), grid, xlabel('t'), ylabel('n_2'), title('dn_2/dt=\alphan_1\cdotn_2-\betan_2'), legend(‘n_{10}=100’, ‘n_{10}=110’, ‘n_{10}=90’)

Как видим, при наблюдается рост популяции, при - уменьшение её численности.

Рис. 5. Зависимости численности особей различных полов от времени в зависимости от начальных условий

Математические модели на основе обыкновенных дифференциальных уравнений 2-ого порядка.

Задача 1: Свободное падение тела

Рассмотрим свободно падающее тело. На него действует лишь сила тяжести. Тогда второй закон Ньютона записывается ,

откуда

.

Для преобразования дифференциального уравнения второго порядка в систему дифференциальных уравнений первого порядка сделаем замену:

, .

Тогда получим

(13)

Алгоритм

Создадим отдельный М-файл, где зададим функцию, описывающую ОДУ

function dy=fall(t,y,g)

dy=[y(2); -g];

Получим решение при помощи оператора решателя ОДУ по схеме Рунге-Кутты 4-5 порядков (для различных значений ускорения свободного падения: на Земле, на Юпитере и на Луне) для тела, падающего с высоты 100 м без начальной скорости в интервале времени от 0 до 15 секунд.

[t,y]=ode45(@fall, [0 15], [100 0],[],9.8);

[t1,y1]=ode45(@fall, [0 15], [100 0],[],24.9);

[t0,y0]=ode45(@fall, [0 15], [100 0],[],1.62);

Графически представим результаты расчетов, разбив графическое окно на два подокна: слева – зависимость координаты от времени, справа – скорости. Обратите внимание на оператор axis, не позволяющий отображать результаты, лежащие ниже нулевого уровня.

subplot(121), plot(t,y(:,1),t1,y1(:,1),’--',t0,y0(:,1),’:’),

grid, xlabel('t'), ylabel('z'), title(‘dz/dt=V_z’)

axis([0 15 0 100])

legend('g=9.8m/s^2’; 'g=24.9 m/s^2’, 'g=1.62 m/s^2’)

subplot(122), plot(t,y(:,2),t1,y1(:,2),’--',t0,y0(:,2),’:’),

grid, xlabel('t'), ylabel('V_z'), title(‘dV_z/dt=-g’)

Результаты представлены на рисунке 6. Как и ожидалось, зависимость высоты от времени представляет собой отрезок параболы, а график скорости – линейную зависимость. Скорость нарастает тем быстрее, чем больше ускорение свободного падения. Проверив себя на данном примере, можем перейти к более сложным задачам.

Рис. 6. Зависимости высоты и скорости тела от времени при свободном падении

Задача 2: Падение тела с учетом вязкого трения

Апробировав алгоритм на простой задаче, аналитическое решение которой нам известно, мы можем сделать шаг «навстречу реальности». Как известно, в атмосфере на падающее тело начинает действовать сила сопротивления воздуха, поэтому

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

, (14)

где - вязкость среды. Тогда уравнение динамики можно переписать в виде

Учитывая, что масса капельки

,

получим:

(15)

Проведя аналогичную предыдущему разделу замену, получим математическую постановку задачи

Алгоритм

function dy=fall1(t,y,g,eta,ro,R);

dy=[y(2);

-g-4.5*eta/ro/R/R*y(2)];

Рис. 6. Зависимости высоты и скорости капель дождя от времени при учёте вязкого трения

Произведем расчет для двух капель воды диаметром 1 мм и 0.2 мм.

R1=5e-4; R2=2e-4;

eta=1.86e-5;

ro=1000;

[t1,y1]=ode45(@fall1,[0 15],[100, 0], [], 9.8, eta, ro, R1)

[t2,y2]=ode45(@fall1,[0 15],[100, 0], [], 9.8, eta, ro, R2)

subplot(121)

plot(t1, y1(:, 1), t2, y2(:,1))

xlabel('Время'), ylabel('Высота'),

axis([0 15 0 100])

subplot(122)

plot(t1, y1(:, 2), t2, y2(:,2))

xlabel('Время'), ylabel('Скорость')

Результаты расчета, показывают, что вязкое трение оказывает слабое действие на большую каплю. Между тем, маленькая капля уже через две секунды прекращает разгон и в дальнейшем двигается равномерно.

Задача 3: Падение тела с учетом турбулентного трения

Условие реализации ламинарного режима обтекания тела определяется числом Рейнольдса , где R – характерный размер, v – скорость, а - вязкость воздуха – оно должно быть менее 100. Учитывая, что вязкость воздуха , ясно, что это возможно только для очень маленьких тел, размером менее 1 мм и двигающихся со скоростью менее 1 м/с – видимо, для маленьких капелек воды и частиц пыли в атмосфере.
Между тем, даже скорость капель дождя составляет 10 м/с, поэтому для прочих тел преобладающим является турбулентное трение

. (16)

Заметим, что коэффициент турбулентного трения в общем случае не может быть определен. Его необходимо определять из условия равновесия силы тяжести и силы сопротивления. Например, для капли дождя диаметром 1 мм

Из уравнений динамики получим выражение для ускорения

(17)

После замены переменных получаем систему уравнений первого порядка

Алгоритм

Создаем функцию в отдельном М-файле

function dy=fall2(t,y,g,gamma);

dy=[y(2);

-g-gamma*(y(2)).^2.*sign(y(2))];

Произведём расчёт, используя эту функцию.

ro=1000;

R=5e-4;

gamma=4*pi*R^3*ro/3*9.8/100;

[t2,y2]=ode45(@fall2,[0 15],[100, 0], [], 9.8,gamma)

subplot(121)

plot(t, y(:, 1), ’-.’, t1, y1(:,1),’:’,t2, y2(:,1))

xlabel('Время t, c'), ylabel('Высота z, м’),

axis([0 15 0 100])

subplot(122)

plot(t, y(:, 2),’-.’, t1, y1(:,2), ’:’, t2, y2(:,2))

xlabel('Время t,c'), ylabel('Скорость V_z, м/с')

legend(‘свободное падение’,’вязкое трение’,’турбулентное трение’)

axis([0 15 -50 0])

На итоговом рисунке представлены результаты расчета для капли диаметром 1 мм по трем моделям – в свободном падении, с учетом вязкого и турбулентного трения. Время падения с высоты 100 м по разным моделям различается более чем вдвое.

Рис. 7. Зависимости высоты и скорости капли диаметром 1 мм от времени в различных приближениях

Двумерные задачи с ОДУ 2-го порядка

Следующим этапом является рассмотрение двумерных задач динамики. Типичным представителем этого класса является баллистическая задача

Баллистическая задача без учёта сопротивления среды

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

и может быть расписано вдоль горизонтальной и вертикальной осей

, (18)

Как видно из уравнений, движение вдоль разных осей никак не связано друг с другом и может быть описано по отдельности, а для решения задачи можно воспользоваться уже имеющимися у нас алгоритмами.

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

Рис. 8. Траектории тел, вылетевших под различными углами с начальной скоростью 100 м/с (без учёта сопротивления воздуха)

g=9.8;

alpha = 45*pi/180; V0=100;

Vx=V0*cos(alpha); Vy=V0*sin(alpha);

T=Vy/g;

[t1,y1]=ode45(@fall,[0 2*T],[0 Vy],[],g);

x1=Vx*t1;

plot(,x1,y1(:,1));

xlabel('x,m');ylabel('y,m');

axis([0 Vx*T 0 g*T*T]);

text(x1(30)y1(30.1),'\alpha = 45^{\circ}')

Баллистическая задача c учётом сопротивления среды

Как было показано в предыдущем разделе, турбулентное трение является основной силой, влияющей на движение тела в атмосфере. Оно описывается формулой (16). Уравнение динамики записывается в виде

Для проекций на оси получим

(19)

где модуль скорости

Как видим из (19), в результате турбулентного трения движение по разным координатам уже нельзя рассматривать независимо, уравнения оказываются связанными т.к. скорость и модуль скорости зависит от обеих компонент. Совершив замену

, , , .

мы получим четырехкомпонентый вектор, для которого уравнения (19) запишутся в виде

Алгоритм

Создаем функцию в отдельном М-файле

function dy = ballist(t,y,g,gamma)

V=sqrt(y(2)*y(2)+y(4)*y(4)));

dy=[y(2); - gamma*v*y(2);y(4);-g-gamma*v*y(4)];

Рис. 9. Траектории тел, вылетевших под различными углами с начальной скоростью 100 м/с с учётом сопротивления воздуха

Дата добавления: 2015-09-29