Непрерывные модели

Систематика моделей
МОДЕЛИ
материальные
идеальные
знаковые
интуитивные
математические
дискретные
непрерывные
стохастические

Характеристика моделей:
Материальные модели – это модели, для работы с которыми необходим эксперимент.
Знаковые модели – это любые модели, использующие знаки.
Интуитивные модели – это модели, опирающиеся на жизненный опыт.
Результатом математического моделирования является теория и выводы из нее.

Процесс моделирования:
1. Определение цели моделирования (что дано и что требуется найти).
2. Определение факторов, которые нужно учесть.
3. Определение точности данных и требуемой точности результата. Сопоставление их.
4. Построение нескольких моделей одного и того же объекта (желательно).
5. Проверка адекватности модели.
*Результат моделирования носит приближенный или гипотетический характер!
** В каждой модели договариваются о предположениях, выдвигают гипотезы.
Адекватность – это соответствие модели реальному объекту. Адекватность проверяется здравым смыслом или сравнением результатов моделирования с практическими результатами.

Задача баллистики.
Задача: из катапульты бросают камень с начальной скоростью v0 под углом α к горизонту.
Требуется определить траекторию полета камня и дальность полета.
Предположения модели:
1. Земля – инерциальная система отсчета.
2. Ускорение свободного падения g постоянно.
3. Кривизной Земли можно пренебречь.
4. Сопротивлением воздуха можно пренебречь.
ТЕМА: ПРИМЕРЫ МОДЕЛЕЙ, ПРИВОДЯЩИХ К
ДИФФЕРЕНЦИАЛЬНЫМ УРАВНЕНИЯМ

x
y
α
v0
vx=v0cosα vy=v0sinα-(gt2)/2
x=tv0cosα (1)
y=tv0sinα-(gt2)/2 (2)

Выражаем t из (1) и подставляем в (2). Получаем уравнение траектории:
Чтобы определить дальность полета, нужно решить уравнение:
х1=0;
формула для определения дальности полета.

Задача баллистики с учетом силы сопротивления воздуха.
Из предположений модели исключаем предположение №4.
Предположения модели:
1. Земля – это инерциальная система отсчета.
2. Ускорение свободного падения g постоянно.
3. Кривизной Земли можно пренебречь.
4. Сопротивлением воздуха можно пренебречь.

- подъемная сила, которая учитывается для асимметричных тел;
- лобовое сопротивление.
Считаем, что камень имеет сферическую форму и поэтому подъемной силой можно пренебречь.
Fл=cSρ(v2/2),
где с – коэффициент лобового сопротивления;
S – площадь поперечного сечения тела;
ρ – плотность воздуха;
v – скорость.
S=πR2, где R – радиус тела.
Fл=c πR2 ρ(v2/2),

Знак «минус» показывает противоположное направление скорости.
i
j
x
y
v
Fл
Fлx
Fлx
vy
vx

Подставляем Fл в уравнение (*):
Разделим обе части уравнения на m:
Перейдем от этого уравнения к скалярным уравнениям:

x(t)≡x
y(t)≡y
vx=x’
vy=y’
ax=x’’
ay=y’’
Чтобы решить эту систему, нужно знать начальные условия.
Начальные условия:
x(0)=0; x’(0)=v0cosα;
y(0)=0; y’(0)=vosinα;
Сменим обозначения:
x→U1;
x’→U2;
y→U3;
y’→U4.

Начальные условия:
U1(0)=0;
U2(0)=v0cosα;
U3(0)=0;
U4(0)=v0sinα.
Решением этой системы будут функции при заданных начальных условиях.

ТЕМА: КРАЕВАЯ ЗАДАЧА
Постановка задачи. Снаряд летит с постоянной скоростью под углом α к горизонту.
Цель находится на расстоянии l.
Угол α мы можем менять.
Вопрос: как стрелять?

Система дифференциальных уравнений:
Граничные условия (начальные условия, которые задаются на концах отрезка):

Дальность полета максимальна при α=45º.
Мы можем решить задачу Коши
при α =45º. Получим некоторую траекторию, которая будет больше, чем нужно.
Уменьшим угол α. Возьмем α1<45º и снова решим задачу Коши.
Получим новую траекторию. Здесь 2 варианта – перелет или недолет.
Маловероятно, что мы попадаем в цель. После этого снова меняем угол α:
x(tk)=l;
|x(tk)-l|≤ε.
При достижении заданной точности ε мы можем прекратить решение задачи
и считать полученное решение окончательным решением краевой задачи.

α=45º
решаем задачу Коши
перелет (траектория больше, чем нужно)
уменьшаем угол α
α1<45º
решаем задачу Коши
новая траектория
увеличиваем угол α
45º>α2>α1
решаем задачу Коши
новая траектория
ДА
НЕТ
Достигнута ли требуемая точность ε?
ОТВЕТ

где:
где
— вектор начальных условий, а
— вектор правых частей.
— вектор, содержащий решение,
Задача Коши в векторной форме:

Задача Коши для системы
нормальных
дифференциальных
уравнений
РЕШЕНИЕ ЗАДАЧИ
?
КАК ЭТО СДЕЛАТЬ В MATHCAD?
Функция rkfixed

Функция rkfixed(
, a, b, n, F) решает задачу Коши методом Рунге—Кутты
с постоянным шагом.
– вектор начальных условий,
F – вектор правых частей уравнений,
n – число узлов (точек разбиения) на отрезке [a,b].
Задача Коши для системы
нормальных
дифференциальных
уравнений
Функция rkfixed
Нижняя граница номеров элементов массивов полагается равной 1
и задается системной константой ORIGIN, которая по умолчанию
равна 0.

Пример. Решить систему дифференциальных уравнений на отрезке [0,3] с шагом 0.1, построить графики найденных функций.
Решение в Mathcad

В этой модели учитывается один показатель: количество публикаций в данный момент.
х(t)≡х;
ПРЕДПОЛОЖЕНИЕ: скорость роста количества публикаций пропорциональна их количеству в данный момент времени.
Модель развития науки.

Запишем дифференциальное уравнение:
Решением уравнения такого вида является функция вида еt; et→∞; t→∞.
! Такую модель можно применять только на конечном участке времени !
Пусть b – уровень насыщения. Тогда дифференциальное уравнение
можно записать следующим образом:

Модель численности популяций.
N(t) – численность популяции в момент времени t.
Скорость v роста популяции прямо пропорциональна ее численности
в данный момент времени.
где μ – разность между коэффициентом смертности и коэффициентом рождаемости.
Если μ>0, то численность популяции неограниченно растет, т.е. N→∞.
! В этой модели не учитываются ограничения среды. !
Пусть b – емкость среды. Тогда дифференциальное уравнение запишется следующим образом:
! В этой модели ограничения среды учтены !

Пример 3. Модель гонки вооружений Ричардсона.
x(t)≡x – расходы на вооружение «зеленых»;
y(t)≡y – расходы на вооружение «желтых».
ПРЕДПОЛОЖЕНИЕ 1: Затраты «зеленых» зависят от уровня затрат «желтых», и наоборот.
Начальные условия:
ПРЕДПОЛОЖЕНИЕ 2: скорость роста затрат на вооружение уменьшается с ростом самих затрат.
Чем больше х, тем меньше dx/dt.

ПРЕДПОЛОЖЕНИЕ 3: Предположим, что у каждого государства могут быть какие-то претензии.
Введем коэффициенты претензий r и s.
r,s≥0 – коэффициенты претензий;
r,s<0 – коэффициент доброй воли.
Система уравнений запишется так:
Начальные условия:

КОГДА ДОСТИГАЕТСЯ РАВНОВЕСИЕ?
Точка равновесия достигается при следующем условии:
x=(a/m)y+(r/m);
x=(n/b)y-s/b;
(a/m)y+r/m= (n/b)y-s/b;
y(a/m-n/b)=-s/b-r/m.

Точка, в которой достигается равновесие, имеет следующие координаты:
Каждая точка этой координатной плоскости характеризуется:
1. Координатами х и у.
2. Моментом времени t.
3. Скоростью изменения координат dx/dt и dy/dt.

- точка равновесия. Стрелка показывает притяжение точек к прямой.
Если dx/dt=0, то прямая G: ay-mx+r=0;
Если dx/dt>0, то x(t) возрастает (точка движется вправо);
Если dy/dt=0, то Z: bx-ny+s=0.
х
y

х
y

Исследуем поведение системы в зависимости от начальных условий.
В зависимости от коэффициентов m, n, r, s
положение точки равновесия может быть различным.
Рассмотрим следующие варианты:
Вариант 1. mn-ab>0; r,s>0.
! При любых начальных условиях есть стремление к точке равновесия!
G: y=(m/a)x-r/a
Z: y=(b/n)x+s/n.
m/a-b/n=(mn-ab)/an>0,
т.е. коэффициент m/a прямой G больше коэффициента прямой Z,
следовательно, угол наклона прямой G больше.

Вариант 2. mn-ab>0; r, s<0.
x(t), y(t)≥0.
Точка равновесия не достижима ни при каких начальных условиях.
Область 1
Область 2
Область 3
G
Z
В областях 1 и 2 к равновесию не приходим.
В области 3 возможно полное разоружение сторон (решение стремится к точке (0;0)).

Вариант 3. mn-ab<0; r,s<0.
В областях 1 и 2 есть стремление к положению равновесия.
В области 3 – бесконечная гонка вооружений.
В области 4 есть стремление к всеобщему разоружению.

ПРОВЕРКА АДЕКВАТНОСТИ МОДЕЛИ делалась на реальных данных.
Рассматривались 2 блока:
х – Россия и Франция;
у – Германия и Австро-Венгрия.
ПРЕДПОЛОЖЕНИЕ: а=b; m=n.
При таком предположении система
перепишется следующим образом:
Складываем два уравнения:
d(x+y)/dt=(a-m)y+(a-m)x+r+s.
Пусть a-m=k, а r+s=f.
z=x+y. - суммарные затраты на вооружение 2-х блоков
dz/dt=kz+f – это уравнение можно использовать
для анализа.

Таблица 1. Затраты на вооружение.

Δz/Δt
z
10
230
250
210
20
30
Вывод: в эту модель хорошо укладываются реальные данные.

Моделирование прогиба балки

На свободный конец балки единичной длины и постоянной жесткости (ЕJ=1, Е – модуль Юнга , J – коэффициент жесткости) действуют вертикальная сила Q и горизонтальная сила P
Здесь S – длина дуги деформированной оси балки от его начала до точки с абсциссой х. Касательная к оси балки в этой точке образует с осью х угол φ(S).

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

Если известно φ(S), то декартовы координаты точек изогнутой оси балки можно определить из следующих соображений:
малому изменению декартовых координат Δх и Δу соответствует малое изменение длины дуги стержня ΔS. Малый отрезок дуги можно заменить отрезком прямой, и тогда Δх= ΔS cosφ, Δу= ΔS sinφ. Отсюда
Задавая различные значения Р и Q, можно построить кривые у(х) при заданных нагрузках.

Решение задачи КОШИ
для обыкновенного дифференциального уравнения
odesolve(x,b,step),
Обращение к функции имеет вид y:=odesolve(x,b) или y:=odesolve(x,b,step),
где y – имя функции, содержащей значения найденного решения,
x — независимая переменная,
b — конец промежутка интегрирования,
step – шаг (можно не указывать).
Перед обращением к функции odesolve надо записать ключевое слово given, затем ввести уравнение и начальные условия. При вводе уравнения и начальных условий используется знак символьного равенства. При записи неизвестной функции в уравнении обязательно требуется указать в скобках аргумент – независимую переменную. Производную можно записывать как , так и (знак производной ставится сочетанием клавиш Ctrl и F7). Аналогично записываются и производные высших порядков.

Решение в Mathcad

Применение метода конечных разностей
дифференциальная задача сводится к решению системы алгебраических уравнений.
Длина балки разбивается на n равных интервалов длиной h.
В каждой точке разбиения записывается исходное дифференциальное уравнение, при этом во всех узлах сетки, кроме нулевого, вторая производная заменяется ее конечно-разностной аппроксимацией:
При i = 0 учет первого краевого условия дает .

При для аппроксимации производной требуется значение угла φn+1 в несуществующем узле с номером n+1. Оно определяется из второго краевого условия:
Отсюда
В итоге имеем систему нелинейных уравнений, решением которой являются значения угла в точках разбиения:

При малых прогибах балки можно считать
и решать систему линейных уравнений относительно
.

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

Общий вид дифференциального уравнения 2–го порядка с двумя независимыми переменными:
где u≡u(x,y), неизвестная функция, х и у – независимые переменные.
Дифференциальные уравнения в частных производных имеют бесконечное множество решений.
Чтобы получить единственное решение для конкретного физического процесса, необходимо задать дополнительные условия: начальные и краевые (граничные). Их вид и количество зависит от типа уравнения.

ТИПЫ уравнений в частных производных
Тип уравнения определяется по виду коэффициентов и соотношениям между ними. Обычно рассматривают уравнения трех основных типов: параболического (D = b2–4ac=0); гиперболического (D>0); эллиптического (D<0).
Если все коэффициенты a, b, c, d, e, f, g не зависят от u, x, y, то получается уравнение с постоянными коэффициентами.
Если g линейно зависит от u, а остальные коэффициенты зависят только от x, y, уравнение называется линейным. Если a = b = c = f, а d ≠ 0, e ≠0, то имеем уравнение 1–го порядка.

Рассмотрим дифференциальное уравнение в частных производных первого порядка:
, u(t,x) – неизвестная функция двух независимых переменных t, x;
u(0,x)=φ(х) – начальное условие при t=0;
u(t,0)=ψ(t) – краевое условие при х=0.
Рассмотрим прямоугольную область: 0≤t≤Т 0≤х≤1. Строим прямоугольную сетку: хj=jh, ti=iτ, i=0,1,2,…,n, j=0,1,…,m. Вместо функций u(t,x), φ(х), ψ(t), f(t,х) будем рассматривать сеточные функции uij, fij, φij, ψij.
Аппроксимируем производные по шаблону на рис.3. a).

Уравнение Пуассона
Это уравнение эллиптического типа, которое можно решить в Mathcad с помощью встроенной функции
Уравнение Пуассона описывает стационарное распределение температуры u(x,t) на плоскости, в которой имеются источники тепла с интенсивностью f(x,y).
Область решения уравнения Пуассона должна иметь форму квадрата, граничные условия задаются на всех сторонах квадрата. В самом простом случае с нулевыми или постоянными граничными условиями (постоянная температура на всей граничной области) можно использовать встроенную функцию multigrid.

Функция multigrid (f, r) имеет два аргумента:
f – имя матрицы, задающей правую часть уравнения; r – параметр численного метода, обычно полагают r = 2. Количество точек разбиения стороны квадрата m должно быть равно степени двойки: 4, 8, 16, и т.д. Графически функцию можно изобразить поверхностью, выбрав на панели инструментов «построение поверхности» и указав в шаблоне имя матрицы, содержащей значения неизвестной функции u(x,y).

Пример 2. Решить уравнение Пуассона, если при , m = 32 , имеются два источника тепла и один сток тепла.
Решение

Аппроксимация частных производных
функция двух переменных.
Пусть область решения заменена сеткой, узлы сетки имеют координаты
(i,j) – целочисленные координаты узла.
При аппроксимации частных производных рассматривают некоторую окрестность узла, называемую шаблоном.
Выбор шаблона обычно определяется конкретной задачей и определенными требованиями к решению. Наиболее популярный шаблон – так называемый «крест» ) и различные его подмножества.

Шаблон «крест»
Шаг по переменной t равен τ.
по переменной x шаг равен h

По этому шаблону частные производные функции
в узле (i,j) можно аппроксимировать следующим образом:
;
;
;
;
;
;
;

Дифференциальное уравнение в частных производных первого порядка:
,u(t,x) – неизвестная функция двух независимых переменных t, x;
u(0,x)=φ(х) – начальное условие при t=0;
u(t,0)=ψ(t) – краевое условие при х=0.
Рассмотрим прямоугольную область: 0≤t≤Т 0≤х≤1. Строим прямоугольную сетку: хj=jh, ti=iτ, i=0,1,2,…,n, j=0,1,…,m. Вместо функций u(t,x), φ(х), ψ(t), f(t,х) будем рассматривать сеточные функции uij, fij, φij, ψij.

Аппроксимируем производные по шаблону на рис
Запишем исходное уравнение в каждом внутреннем узле области решения:

Отсюда
где i=0,1,2,…,n; j=0,1,2,…, m
Вычисления идут по слоям. На нулевом слое u0j=φj (φi вычисляется как φ(xj)), ui0=ψ(ti).
Получилась явная разностная схема: значения сеточной функции в каждом узле верхнего слоя вычисляются через значения на предыдущем слое.

Схема устойчива при 0≤τ≤h/a, при этом u(t,x), φ(х), ψ(t) должны быть дважды непрерывно дифференцируемы; f(х,t) имеет непрерывные первые производные. Сеточное решение сходится к точному при h→0; τ→0. При а<0 условие устойчивости не выполняется, и схема не сходится.
При а<0 сходящаяся разностная схема получается при использовании шаблона, показанного на рис.
Эта схема устойчива при а<0, если выполнено условие τ≤–h/a. При а>0 эта схема не сходится.

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

Пример.

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

x
u