Задачи: краевая, о собственных числах, краевая для уравнения Пуассона, Коши, теплопроводности и линейного программирования
Выбери формат для чтения
Загружаем конспект в формате pdf
Это займет всего пару минут! А пока ты можешь прочитать работу в формате Word 👇
Оглавление
Лекция 1 ............................................................................................................................................................................. 4
КРАЕВАЯ ЗАДАЧА ......................................................................................................................................................... 4
Практическая работа 1 .................................................................................................................................................... 19
Лекция 2 ........................................................................................................................................................................... 20
ЗАДАЧА О СОБСТВЕННЫХ ЧИСЛАХ ...................................................................................................................... 20
Практическая работа 2 .................................................................................................................................................... 30
Лекция 3 ........................................................................................................................................................................... 31
3. КРАЕВАЯ ЗАДАЧА ДЛЯ УРАВНЕНИЯ ПУАССОНА .......................................................................................... 31
Практическая работа 3 .................................................................................................................................................... 39
Лекция 4 ........................................................................................................................................................................... 40
ЗАДАЧА КОШИ (задача с начальными условиями) ................................................................................................... 40
Практическая работа 4 .................................................................................................................................................... 63
Лекция 5 ........................................................................................................................................................................... 64
ЗАДАЧА ТЕПЛОПРОВОДНОСТИ .............................................................................................................................. 64
Практическая работа 5 .................................................................................................................................................... 73
Лекция 6 ........................................................................................................................................................................... 74
ЗАДАЧИ ЛИНЕЙНОГО ПРОГРАММИРОВАНИЯ .................................................................................................... 74
Практическая работа 6 .................................................................................................................................................... 82
Лекция 1
КРАЕВАЯ ЗАДАЧА
Введение
1.1. Обыкновенные дифференциальные уравнения
1.2. Общий вид краевой задачи на отрезке (a,b) для
дифференциального уравнения второго порядка
1.3. Метод конечных разностей
1.4. Порядок аппроксимации производных
1.5. Решение краевой задачи в MATLAB
Практическая работа 1
Введение
В
настоящее
время
применение
численных
методов
и
методов
компьютерного моделирования не ограничивается только естественнонаучными
областями знания, но и активно востребовано в сфере, традиционно относимой
к гуманитарной. Это - история, политика, экономика, психология и некоторые
другие области знаковой деятельности человека.
Активное продвижение математических методов и технологий во все сферы
знаковой деятельности человека стало возможным с появлением компьютера,
т.е. с появлением компьютерной или, более точно, информационной индустрии.
Традиционное использование в науке компьютера, как универсального
вычислительного устройства получило в последнее время новое содержание в
рамах таких понятий, как виртуальная реальность и интерактивный интерфейс.
Современные вычислительные мощности и средства программирования
позволяют в ряде случаев не просто решить ту или иную задачу, но и построить
нечто большее - виртуальную реальность, в которой исследователь может
осознать решение задачи во всем возможном многообразии. Привлечение
средств
интерактивного интерфейса позволяет активно манипулировать
доступным многообразием для получения оптимального в том или ином смысле
решения исходной задачи. В этом контексте курс лекций построен таким
образом, чтобы студенты могли непосредственно использовать такой широко
известный пакет программирования, как MATLAB, для активного овладения
вычислительными методами.
В курсе строительной информатики будут рассмотрены современные
методы численного и компьютерного моделирования для решения прикладных
задач строительной отрасли и их реализациями в системе MATLAB.
1.1 Обыкновенные дифференциальные уравнения
Обыкновенные дифференциальные уравнения занимают видное место в
науке. Они широко используются в механике при расчете различного рода
движений, в других разделах физики, например, при расчетах электрических
цепей, в химической кинетике при анализе баланса реагентов, в популяционной
биологии (модели типа хищник-жертва), в машиностроении, в строительстве при
расчете сооружений, в экономике, в политике и во множестве других областях.
В общем и целом обыкновенные дифференциальные уравнения являются
одним из самых востребованных инструментов описания и расчета самых
разнообразных явлений и процессов в науке, технике и обществоведении.
Обыкновенное
дифференциальное
уравнение
-
это
уравнение,
содержащее аргумент, искомую функцию этого аргумента и ее производные
различных порядков:
F ( x, y, y' , y' ' ,.., y ( n ) ) 0
(1.1)
Порядок дифференциального уравнения — это наивысший порядок
производной искомой функции в данном уравнении.
Решением дифференциального уравнения называется функция, которая
при подстановке в дифференциальное уравнение обращает последнее в верное
равенство.
Частным случаем (1.1) является дифференциальное уравнение первого
порядка, разрешенное относительно первой производной:
y'=f(x,y)
(1.2)
Если неизвестная функция y является вектором, то тогда задача сводится к
системе дифференциальных уравнений.
При этом известно, что дифференциальные уравнения и системы
произвольного порядка могут быть сведены к системам дифференциальных
уравнений первого порядка, используя введение вспомогательных функций,
количество которых соответствует порядку дифференциального уравнения.
Например, для дифференциального уравнения второго порядка
y''=f(x,y,y')
можем ввести новые функции:
y1 y,
тогда
дифференциальное
y2 y1' y ' ,
уравнение
второго
порядка
сводится
к
дифференциальным уравнениям первого порядка
y1 ' y2
y2 ' f ( x, y1 , y2 ).
Известно, что дифференциальное уравнение n-го порядка имеет множество
решений, которые в общем случае зависят от n параметров. Для выделения
единственного решения необходимо наложить n дополнительных условий.
1.2. Общий вид краевой задачи на отрезке (a,b) для дифференциального
уравнения второго порядка
y P x y q x y f x , x a, b ,
1 y a 1 y a 1
y b y b краевые условия.
2
2
2
(1.3)
Здесь все величины, кроме y x , предполагаются заданными.
Отметим еще раз, что существует некоторое бесконечное множество
функций, удовлетворяющих приведенному дифференциальному уравнению, и
лишь наличие области изменения x ( x (a, b) – геометрия конструкции) и
условий на краях области обеспечивает единственность решения. Однако и это
решение (хотя оно и единственное) в общем случае нельзя представить в
аналитической форме. Поэтому, как правило, задача решается численным
методом.
1.3. Метод конечных разностей
Это наиболее простой и при этом достаточно эффективный способ решения
краевой задачи. Суть его состоит в следующем.
Разобьем отрезок (a, b) на ( N 1) -отрезков (рис.1.3).
Рис. 1.3.
Введем обозначения:
x i – координаты точек разбиения,
i – номер точки разбиения ( i 1, 2, ... , N ),
hi xi 1 xi – длина i -го отрезка (шаг разбиения),
x xi 1
~ h h
hi i i 1 i 1
– «средний» шаг,
2
2
y i y xi , при этом y1 y a ,
y n y b ,
Pi P xi ,
qi q xi ,
f i f xi .
Производные в i -ой точке заменим разностными соотношениями:
yi
y i
yi
yi 1 yi
– правая разность,
hi
y i y i 1
– левая разность,
hi 1
(1.4)
yi 1 yi 1
– центральная разность.
~
2hi
При h 0 , в соответствии с определением производной, все три величины
будут стремиться к y ( xi ) . Использование той или другой из них зависит от
конкретной ситуации.
Вторая производная в i -ой точке может быть приближенно представлена
разностным отношением в виде разности от первых, например,
y y i y i y i 1 ~
hi ,
y i i 1
hi
hi 1
при hi h const
y i
y i 1 2 y i y i 1
h2
(1.5)
Здесь предложена наиболее часто употребляемая формула второй разности.
Пользуясь приведенными обозначениями и формулами, можно представить
задачу (1.3) в каждой i -ой точке относительно величин y i следующим образом:
y 2 y1
y
1 , i 1,
1
1
1
h1
y i 1 y i y i y i 1 ~
y y
hi Pi i 1 ~ i 1 qi y i f i , i 2,3,..., N 1 (1.8)
hi 1
2hi
hi
y y N 1
2, i N.
2 y N 2 N
h N 1
Это и есть разностный аналог краевой задачи (1.3). Здесь N неизвестных y i
и N уравнений. Приводя подобные члены, получим:
a11 y1 a12 y 2 b 1
i 1,
..............................................
ai , i 1 y i 1 ai ,i yi ai , i 1 y i 1 b i , i 2,3,..., N 1,
..............................................
a N , N 1 y N 1 a N , N y N bN
i N.
(1.6)
или в матричной записи
Ay b ,
(1.7)
где
a11
a
21
0
A
0
a12
a 22
a 23
a i ,i 1
a i ,i
a i ,i 1
y1
b1
y
b
2
2
, y , b ,
yi
bi
a N 1, N
a N , N
yN
b N
a N , N 1
при этом элементы матрицы и вектора правой части определяются по
следующим формулам:
a11 1
1
p
a i ,i 1
i
hi 1 2
i 2,3,..., N 1,
1
h1
, a12
1
h1
,
b1 1 ,
~
1 p ~
2
/ hi , a ii qi
, a i ,i 1 i / hi , bi f i ,
2
hi hi 1
hi
a N , N 1
2
hN 1
, a N ,N 2
2
hN 1
, bN 2 .
Отметим, что A – трехдиагональная матрица.
После решения системы (1.8) получим приближенное решение задачи (1.3)
в i -х точках. Соединяя y i ломаной линией, получаем приближенное решение
y x во всех точках области (a, b) .
Точность решения задачи (1.3) зависит от точности аппроксимации
значений производных разностями, которые, в свою очередь, зависят от
величины шага разбиения.
Рассмотрим краевую задачу.
40
40
y 1 4 x (1 x ) y 2.4( 2 1 4 x (1 x ) x (1 x ));
y (0) 0
краевые условия
y (1) 0
x (0,1),
Решение на MATLAB
матрица системы разностных уравнений
1.000
0.000
0.000
0.000
0.000
1.000 -2.435
1.000
0.000
0.000
0.000
1.000 -2.357
1.000
0.000
0.000
0.000
1.000 -2.323
1.000
0.000
0.000
0.000
1.000 -2.313
0.000
0.000
0.000
0.000
1.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
1.000
-2.323
1.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
1.000
-2.357
1.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
1.000
-2.435
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
1.000
1.000
вектор правой части системы разностных уравнений
0.000
-0.189
-0.236
-0.256
-0.263
-0.256
-0.236
-0.189
0.000
0.450
0.562
0.600
0.562
0.450
0.262
0.000
решение y:
0.000
0.262
График решения на рисунке 1.4
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Рис. 1.4 График решения краевой задачи
1.3. Порядок аппроксимации производных
Для упрощения оценок предположим, что hi=h=const. По формуле
Тейлора:
f i 1
h2
h3
h 4 4
f xi h f xi hf xi
f xi
f xi
f xi O h 5 ,
2
6
24
где
Ch , C const .
O hk
k
Тогда для правой разности получим
f i 1 f i
h
f x i
f x i hf x i
h2
f x i ... f x i
2
h
h
f x i ... f x i O h аппроксимация первого порядка.
2
Аналогично для левой разности
f i f i 1
h
f x i f x i ... f x i O h ,
h
2
т.е. также имеет место первый порядок аппроксимации.
Рассмотрим также центральную разность
f i 1 f i 1
2h
f x i hf x i
f xi hf xi
h2
h3
f xi
f x i ...
2
6
2h
h2
h3
f xi
f xi ...
2
6
2h
h2
f xi
f xi ... f xi O h 2
6
–
аппроксимация
второго
порядка.
Следовательно, центральная разность, как правило, является более точной
аппроксимацией первой производной.
Установим порядок аппроксимации для второй производной
f i 1 2 f i f i 1
h2
h2
h3
h 4 4
f xi hf xi
f xi
f xi
f xi ...
2
6
24
h2
f x
2 2i
h
2
h
h3
h 4 4
f xi hf xi
f xi
f xi
f xi ...
2
6
24
h2
h 2 4
f x i
f x i ... f x i O h 2 – аппроксимация второго порядка.
12
Подобным образом можно установить порядок аппроксимации для любой
разностной формулы.
1.4. Решение краевой задачи в MATLAB
В системе MATLAB существуют функции для решения дифференциальных
уравнений.
Для решения обыкновенных дифференциальных уравнений первого
порядка (1.11) в системе MATLAB используется функция bvp4c.
y ' f ( x, y ) ,
(1.11)
где
y-
искомая вектор-функция;
y' -
покомпонентная производная вектор-функции y по переменной х;
f ( x, y ) - заданная вектор-функция правых частей;
x [ a , b] .
Решение задачи ищется в форме сеточной функции: отрезок [a,b] делится
точками х1=a>A = gallery('lehmer',4)
A =
1.0000
0.5000
0.3333
0.2500
0.5000
1.0000
0.6667
0.5000
0.3333
0.6667
1.0000
0.7500
0.2500
0.5000
0.7500
1.0000
Вычислим собственные значения матрицы A. Результат будет в виде векторстолбца.
>>e = eig(A)
e =
0.2078
0.4078
0.8482
2.5362
Результат можно получить и в виде диагональной матрицы, указав
параметр 'matrix'.
>>D = eig(A,'matrix')
D =
0.2078
0.4078
0.8482
2.5362
Получить собственные значения и собственные вектора можно, задав
два выходных параметра.
>> [x,D]=eig(A)
x =
0.0693
-0.3618
0.7694
-0.5219
-0.4422
0.7420
0.0486
-0.5014
-0.8105
-0.1877
0.3010
0.4662
0.3778
0.5322
0.5614
0.5088
0.4078
0.8482
2.5362
D =
0.2078
2.4 Алгоритм определения собственных значений оператора краевой
задачи на MATLAB
Рассмотрим алгоритм решения задачи для определения минимального
собственного значения оператора краевой задачи
Ry Py, x (0,1) ,
y (0) 0
y (1) 0 краевые условия,
где R(x)=7.5x(1-x).
Система конечно-разностных уравнений этой краевой задачи в матричном
виде:
Ay PBy .
Алгоритм решения задачи определяется последовательностью действий:
1. Выделить память под матрицы с одновременным их обнулением через
функцию zeros.
2. Сформировать трехдиагональную матрицу А и диагональную
матрицу B.
3. Получить собственные значения и собственные вектора для B 1 A
через функцию eig.
4. Вывести полученные результаты.
Решение задачи состоит в составлении системы при разбиении отрезка [ 0 , 1
] на ( N 1 )-частей с постоянным шагом h 1 / ( N 1)
7.5x(1 x ) y Py, x (0,1) ,
y (0) 0
y (1) 0 краевые условия.
Результат работы алгоритма на MATLAB
Введите число точек разбиения n= 7
Матрица А
2.00 -1.00 0.00 0.00 0.00 0.00 0.00
-1.00 2.00 -1.00 0.00 0.00 0.00 0.00
0.00 -1.00 2.00 -1.00 0.00 0.00 0.00
0.00 0.00 -1.00 2.00 -1.00 0.00 0.00
0.00 0.00 0.00 -1.00 2.00 -1.00 0.00
0.00 0.00 0.00 0.00 -1.00 2.00 -1.00
0.00 0.00 0.00 0.00 0.00 -1.00 2.00
Матрица B
0.02 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.01 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.01 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.01 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.01 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.01 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.02
Матрица форм потери устойчивости Y
0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.03 0.11 -0.24 -0.38 -0.47 -0.21 -0.42
-0.20 -0.44 0.54 0.33 -0.13 -0.36 -0.48
0.49 0.55 -0.03 0.49 0.34 -0.45 -0.30
-0.66 -0.00 -0.54 -0.00 0.54 -0.48 0.00
0.49 -0.55 -0.03 -0.49 0.34 -0.45 0.30
-0.20 0.44 0.54 -0.33 -0.13 -0.36 0.48
0.03 -0.11 -0.24 0.38 -0.47 -0.21 0.42
0.00 0.00 0.00 0.00 0.00 0.00 0.00
Вектор критических сил P
420.0000
315.0000
225.0000
150.0000
90.0000
15.0000
45.0000
Для отображения графически форм потери устойчивости добавим в М-файл
код
x=0:h:1;
for i=1:5
subplot(7,1,i);
plot(x,Y(:,i));
end
subplot(7,1,6);
plot(x,Y(:,7))
subplot(7,1,7);
plot(x,Y(:,6))
Результат представлен на рисунке 2.3. Из рисунка видно, что для некоторых
форм для отображения в графическом виде недостаточно того количества точек,
что использовались при решении задачи.
Рис. 2.3 Формы потери устойчивости
Для гладкости отображения
можно воспользоваться интерполяцией
кубическими сплайнами. Сплайн - это непрерывная гладкая функция, которая
на отрезках области определения равна полиномам определенной степени, т.е.
при таком способе интерполяции точки попарно соединяются отрезками
полиномов. На практике используют полиномы невысокой степени, чаще —
третьей, т.е. имеет место интерполяция кубическими сплайнами. Таким образом,
аппроксимация происходит не глобально на всем интервале, а по отдельности на
каждом частичном интервале между соседними узлами. График сплайна
проходит через узлы, в которых обеспечивается стыковка - совпадение значений,
совпадение производных для отсутствия изломов, иногда обеспечивается и
совпадение вторых производных и т.д. Могут задаваться и краевые условия на
концах интервала.
Для
выполнения
интерполяции
кубическими
сплайнами
можно
использовать функцию spline, синтаксис которой имеет вид:
yy=spline(x,y,xx)
В результате будет выполнена интерполяция значений вектора y, заданного
для значений аргумента, представленного вектором x. Функция возвращает
вектор yy, содержащий значения интерполирующей функции для значений
аргумента, заданного вектором xx.
Для нашей задачи интерполяция кубическими сплайнами для значений,
заданными вектором x
и
столбцом матрицы
Y(:,i)
может быть
представлена следующим фрагментом М-файла:
x=0:h:1;
h=h/10;
xx=0:h:1;
for i=1:5
subplot(7,1,i);
yy=spline(x,Y(:,i),xx);
plot(xx,yy);
end
subplot(7,1,6);
yy=spline(x,Y(:,7),xx);
plot(xx,yy)
subplot(7,1,7);
yy=spline(x,Y(:,6),xx);
plot(xx,yy)
Полученный результат представлен на рисунке 2.4.
Формы потери устойчивости выведены в соответствии со значениями
критических сил: от самой максимальной до самой минимальной сверху вниз.
Рис. 2.4 Формы потери устойчивости.
Практическая работа 2
Лекция 3
3. КРАЕВАЯ ЗАДАЧА ДЛЯ УРАВНЕНИЯ ПУАССОНА
Введение
3.1. Общий вид формулировки краевой задачи для уравнения Пуассона
3.2. Решение задачи Дирихле методом конечных разностей
3.3. Решение задачи Дирихле методом конечных разностей в MATLAB
Практическая работа 2
Введение
Краевая задача для уравнения Пуассона является математическим
описанием разнообразных технических задач: задачи изгиба мембраны, задачи
стационарного распределения температуры в конструкциях, задачи кручения
стержня, входит главной составной частью в задачи теории упругости и т.д.
3.1. Общий вид формулировки краевой задачи для уравнения
Пуассона
Пусть краевая задача рассматривается в относительно некоторой области
(рис. 3.2).
Рис. 3.2
где
( 1 , 2 ) – вектор нормали к границе области :
1 cos( , x ) , 2 cos( , y ) ,
– граница области .
В зависимости от условий на краях области (краевых условий) краевая
задача может иметь разные названия. В частности:
2U F ,
U g
x, y
2U F ,
U
v f
x, y
– задача Дирихле (первая краевая задача),
– задача Неймана (вторая краевая задача),
где
U
U
U
v1
v2
– производная по нормали.
v
x
y
Если на одной части границы заданы условия задачи Дирихле, а на другой
– условия задачи Неймана, тогда таким образом сформулированная задача
называется смешанной краевой задачей.
3.2. Решение задачи Дирихле методом конечных разностей
Рассмотрим задачу Дирихле в прямоугольной области (рис. 3.3).
Рис. 3.3
Тогда ее математическую формулировку можно представить в виде:
2U F ,
0 x 1 , 0 y 2 ,
1 y , x 0
y , x
2
1
краевые
U ( x, y ) x , y 0
3
4 x , y 2 ,
условия.
(3.1)
Здесь все величины заданы, кроме функции U(x,y) – искомая функция.
Примечание. При решении задачи на ЭВМ естественным является
направление оси y вниз, поскольку это соответствует естественному порядку
вывода (печати) на экран и принтер (вывод информации последовательно сверху
вниз).
Разобьем исходную прямоугольную область на мелкие прямоугольники с
помощью сетки с шагом h1 по 1-му направлению (по координате x ) и с шагом
h2 – по 2-му направлению (по координате y ):
h1
1
N1
и h2
2
N2
где N 1 1 и N 2 1 – число узлов сетки по 1-му и 2-му направлению,
соответственно (рис. 3.4).
Будем рассматривать задачу (3.1) в узлах сетки. Каждый узел сетки имеет
номер, определяемый двумя величинами (i, j ) :
j 0,1, ... , N 1 – номер узла по направлению оси x ,
i 0,1, ..., N 2 – номер узла по направлению оси y
Рис. 3.4
Обозначим
U ij U ( x j , yi ) , Fij F ( x j , yi ) , ij ( x j , yi ) .
Вторые производные для внутренних узлов сетки (i, j ) по каждому
направлению заменим вторыми разностями, т.е.
2U ( x j , yi ) U i , j 1 2U ij U i , j 1
,
x2
h12
2U ( x j , yi ) U i 1, j 2U ij U i 1, j
.
y2
h22
Тогда краевая задача в конечных разностях примет вид
Для внутренних узлов
U i , j 1 2U ij U i , j 1
h12
U i 1, j 2U ij U i 1, j
h22
Fij
(3.2)
j 1,2,..., N 1 1 , i 1,2,...N 2 1 ,
1 y i , j 0 ,
y , j N ,
1
2 i
для граничных узлов: U ij ij
x , i 0 ,
3 j
4 x j , i N 2 .
Т.е. для определения U ij во внутренних точках достаточно решить систему
из ( N 1 1) ( N 2 1) уравнений c ( N 1 1) ( N 2 1) неизвестными, используя
заданные значения U ij на границе.
С алгоритмической точки зрения наиболее простыми способами решения
такой системы являются итерационные, в частности, метод Зейделя и метод
простой итерации. Для этого, выделив из уравнений (3.2) диагональный член U ij
, получим
U i , j 1 U i , j 1 U i 1, j U i 1, j
2
2
.
U ij
F
/
ij 2
2
2
2
h
h
h
h
1
1
2
2
Далее, задав начальное приближение U ij0 0 для внутренних точек и
присвоив заданные значения U ij ij для граничных точек, реализуем алгоритм
метода Зейделя в виде
U i,kj1 U i,kj11 U ik1, j U ik1,1j
2
2
U ij
F
/
ij 2
2
2
2
h1
h2
h1 h2
j 1,..., N 1 1 , i 1,..., N 2 1 , k 1, 2, ...
k
(3.3)
Здесь и в дальнейшем k 1, 2, ... – номер итерации.
Счет ведется до тех пор, пока не будет выполнено условие
zk
N 2 1 N1 1
U ijk U ijk 1
i 1
,
j 1
(3.4)
где
– заданная точность,
z k – оценка погрешности на k -ом шаге итерации.
Примечание. Сходимость метода Зейделя для задач такого типа
доказывается в более подробных курсах вычислительной математики. Она
следует из, так называемой, положительной определенности оператора
(матрицы) краевой задачи. Отметим, что в данном случае диагонального
преобладания нет (и не требуется).
3.3. Решение задачи Дирихле методом конечных разностей в MATLAB
Задача
Решить задачу Дирихле для уравнения Пуассона в прямоугольной
области (рис.3.3) используя метод конечных разностей
Для области, имеющей в плане форму квадрата со стороной dl1=dl2=1,
найти форму кровли U(x,y) в двух вариантах:
• подвесная с удельной нагрузкой с=8(g+s);
• надувная с избыточным давлением - с.
Сверху противоположные стороны попарно ограничены кривыми:
U= 4*g/dl2*(dl2-y)*y
U=4*s/dl1*(dl1-x)*x
Принять модель U c знак плюс для подвесной и минус для надувной
кровли соответственно.
Краевыми условиями служат верхние границы стен.
Для g=3 и s=12
Исходные данные
F ( x, y) c , где c 120 ;
1 ( y ) 2 ( y ) 12 y ( 2 y ) , ( x 0 x 1 ) 0 y 1
U ( x, y )
;
3 ( x ) 4 ( x ) 48 x( 1 x ) , 0 x 1 ( y 0 y 2 )
1 2 1;
Результаты приведены ниже, изображения на рисунках 3.5 и 3.6
Введите n1= 8
Введите n2= 6
Введите dl1= 1
Введите dl2= 1
Количество итераций k= 40
Решение u
0.00 5.25 9.00 11.25 12.00
1.67 1.73 2.46 3.11 3.36
2.67 0.36 -0.63 -0.97 -1.05
3.00 0.00 -1.53 -2.20 -2.39
2.67 0.36 -0.63 -0.97 -1.05
1.67 1.73 2.46 3.11 3.36
0.00 5.25 9.00 11.25 12.00
11.25 9.00
3.11 2.46
-0.97 -0.63
-2.20 -1.53
-0.97 -0.63
3.11 2.46
11.25 9.00
5.25
1.73
0.36
0.00
0.36
1.73
5.25
0.00
1.67
2.67
3.00
2.67
1.67
0.00
Рис. 3.5 Подвесная кровля
Рис.3.6 Надувная кровля
Практическая работа 3
Лекция 4
ЗАДАЧА КОШИ (задача с начальными условиями)
Введение
4.1. Метод Эйлера
4.2 Задача об изгибе консоли
4.3 Оценка погрешности метода Эйлера в MATLAB
4.4 Понятие об устойчивости разностной схемы
4.5 Задача Коши для обыкновенных дифференциальных уравнений.
Динамическая система
4.6 Средства MATLAB для решения задачи Коши для обыкновенных
дифференциальных уравнений
4.7 Решение задачи изгиба консоли функцией MATLAB
Практическая работа 4
Введение
Технические задачи, связанные с дифференциальными уравнениями,
включают в себя, как правило, описание области (геометрии конструкции,
временной промежуток и т.д.) При этом все дополнительные граничные условия
могут задаваться только в начале области (начальный момент времени или
основание конструкции). Такая задача в отличие от краевой задачи называется
задачей Коши (задачей с начальными условиями). В частности, это задачи,
связанные с временными процессами (распределение температуры, колебания
конструкции и т.д.).
Пример 1. Пусть задано дифференциальное уравнение 1-го порядка
y F ( x, y )
(4.1)
Будем рассматривать его в области x 0 . Функция F ( x, y ) задана.
Неизвестны функции y(x) и y (x ) .
Как правило, любое дифференциальное уравнение имеет множество
решений. Поэтому для единственности решения требуется дополнительное
условие, например, вида (рис. 4.1):
y (0) y 0 g – начальное условие (задано).
Таким образом, поставленная задача называется задачей Коши и имеет
единственное решение.
Рис.4.1
Пример 2. Приведем пример задачи Коши для уравнения 2-го порядка:
Ry cy f , x 0 ,
y (0) y0 g1
y (0) y g начальные
2
условия.
(4.2)
Неизвестной является функция y(x) .
То есть, для функции y(x) заданы начальное значение y 0 и начальный угол
наклона y 0 (рис. 4.2).
Рис. 4.2
Сведение задачи (4.2) к системе первого порядка
Введем дополнительную функцию
z( x ) y ( x ) .
Подставляя ее в (4.2) получим:
Rz cy f
y z 0 , x 0
y (0) y 0 g1 начальные условия ,
z (0) z0 g 2
(4.3)
где z 0 y 0 .
4.1. Метод Эйлера
В общем случае задача Коши аналитического решения не имеет, поэтому
применяются численные методы, простейшим из которых является метод
Эйлера. Воспользуемся для решения задачи Коши методом конечных разностей.
Для этого зададимся достаточно малым шагом h . Будем рассматривать решение
в точках с координатами xi :xi i h , i 0,1, 2, ... .
Рис.4.3
hi=xi+1-xi
Обозначим y i y ( xi ) ,
f i f ( xi ) ,
y i y ( xi )
и т.д.
Заменим производную конечно-разностным отношением: yi
yi 1 yi
.
h
Получим конечно-разностные уравнения вида:
для примера 1:
y0 g ,
yi 1 yi
F ( xi , yi ) , i 0,1,2,...
h
или
y0 g ,
y i 1 y i h F ( xi , y i ),
i 0,1,2, .
Аналогично для примера 2 (4.3) получим:
y 0 g1
– начальные значения
z0 g 2
z i 1 z i
cy i f i
Ri
h
, i 0,1, 2, ...
y
y
i
1
i
zi
h
или
(4.4)
y 0 g1
– начальные значения
z0 g 2
z i 1 z i h ( cy i f i ) / Ri
, i 0,1, 2, ...
y i 1 y i h z i
(4.5)
Вместо последней формулы можно записать более точную
yi 1 yi h
zi zi 1
.
2
Таким образом, в каждой точке разбиения количество неизвестных
совпадает с количеством уравнений. Приведенные формулы (4.4) и (4.5)
являются рекуррентными формулами вычислений: зная значения xi , yi , zi ,
вычисляем по ним следующие значения xi 1 , yi 1 , zi 1 и т.д. Пользуясь таким
алгоритмом, последовательно получим решение для любого количества точек
разбиения. Вышеизложенный конечно-разностный подход для решения задачи
Коши называется методом Эйлера.
Примечание. В сложных практических случаях применяются различные
модифицированные алгоритмы, связанные с уточнением шага разбиения на
каждом
шаге
пересчета.
Среди
такого
рода
модификаций
наиболее
употребляемыми являются методы типа Рунге-Кутта, излагаемые в специальных
пособиях.
4.2 Задача об изгибе консоли
Рассмотрим задачу об изгибе консоли, закрепленной в левой части.
Y(x)
Рис. 4.4
нелинейное дифференциальное уравнение прогиба оси консоли и
начальные условия
M ( x)
y
(
x
)
[1 ( y ( x)) 2 ]3 ,
EJ ( x)
y (0) y (0) 0 начальные
(4.6)
условия ,
где M (x) – момент, возникающий в консоли от приложенных нагрузок
M ( x ) m P (1 x ) ( x ) q( )d .
x
В частности, при постоянной нагрузке q(x ) : q( x) q const , выражение
для M (x) может быть записано проще (в этом случае легко вычисляется
интеграл в правой части)
( x ) 2
M ( x ) m P( x ) q
2
Обозначив z ( x ) y ( x ) , перейдем к системе двух уравнений 1-го порядка
M ( x)
[1 z 2 ( x )]3
z ( x)
EJ ( x )
,
y ( x ) z ( x )
z (0) 0
– начальные условия.
y (0) 0
Зададим шаг разбиения h и будем рассматривать решение в точках с
координатами xi : xi i h , i 0,1, 2, ... .
уравнения
z0 0 ,
y 0,
0
zi 1 zi M i
(1 zi2 ) 3 , i 0,1, 2, ... .
EJ i
h
y y
i
zi ,
i 1
h
Получим следующие разностные
Этим уравнениям соответствуют рекуррентные формулы пересчета по
методу Эйлера
z0 y0 0
Mi
2 3
zi 1 zi h EJ (1 zi )
i
h
yi 1 yi ( zi zi 1 )
2
,
i 0,1, 2, ...,N 1
Для расчета возьмем следующие соотношения:
M( x )
1
[ 1 ( cx ) ]
2
3
,
Текст М-функции приведен ниже.
Результаты приведены ниже
введите
введите
введите
введите
введите
dl 1
N 10
x0 0
y0 0
z0 0
dl =
1
N =
10
x0 =
y0 =
z0 =
h =
0.1000
x
y
z
EJ c 1, c 0,3, l 1, .
0.1000
0.2000
0.3000
0.4000
0.5000
0.6000
0.7000
0.8000
0.9000
1.0000
0.0015
0.0060
0.0135
0.0240
0.0375
0.0540
0.0735
0.0960
0.1215
0.1500
0.0300
0.0600
0.0900
0.1200
0.1500
0.1800
0.2100
0.2400
0.2700
0.3000
4.3 Оценка погрешности метода Эйлера в MATLAB
В качестве примера рассмотрим поведение оценки погрешности на примере
численного решения уравнения вида:
u ' xu, u(0) 1,
x [0,1].
Это уравнение имеет точное решение
x2
2
ue ,
которое сравним с приближенным решением yn, определяемым на
равномерной сетке xn (n 1) / h, n 1,2,..., n, где h 1/(N - 1) .
Изучим зависимость величины константы в точке x=1
M (1)
1
| y N e1 / 2 |
h
от шага сетки.
Ниже приведен текст М-файла
%Программа изучения погрешности численного решения
%уравнения u'=f(x,u) методом Эйлера
clear
%Задаем набор сеток
Mesh=10:10:10000;
%Организуем цикл расчета методом Эйлера на разных сетках
for k=1:length(Mesh)
%Определяем параметры сетки
N=Mesh(k); h=1.0/(N-1);
%Задаем начальные условия
y(1)=1;
%Алгоритм метода Эйлера
for n=1:(N-1)
y(n+1)=y(n)+(n-1)*h^2*y(n);
end
%Вычисляем М(1) в оценке погрешности численного решения
M(k)=abs(y(N)-exp(0.5))/h;
step(k)=h;
end
Рисуем зависимость М(1) от шага разностной сетки
semilogx(step,M);
На рисунке 4.5 видно, что при стремлении шага сетки к нулю величина М(1)
ведет себя как константа, т.е. не зависит от шага сетки.
Рис.4.5 Изучение погрешности метода Эйлера.
4.4 Понятие об устойчивости разностной схемы
Рассмотрим простейшую задачу Коши вида
y y , x 0
, где 0 .
y ( 0) y 0
В этом случае известно аналитическое решение
y y 0 e x ,
из которого видно, что
y( x)
0 .
x
Разностная формула метода Эйлера имеет вид
yi 1 yi hyi (1 h ) yi .
Откуда следует, что
y i (1 h ) i y 0 .
Обозначим
q 1 h
Тогда
yi q i y 0
При этом, очевидно, что
yi i
0
только при q 1 , что соответствует соотношениям
1 1 h 1 или 0 h 2
откуда следуют ограничения для величины h :
0h
2
.
Последнее неравенство называется условием устойчивости счета.
Абсолютно устойчивая схема
Существуют разностные схемы, обеспечивающие устойчивость счета при
любом шаге h . В частности, для нашего примера абсолютно устойчивой будет
схема
yi 1 yi
y yi
i 1
,
h
2
или
yi 1
1
1
h
2 y
i
h
2
В этом случае
q
1
h
2 .
h
1
2
Очевидно, что
q 1 и yi i
0 ,
т.е. такая схема устойчива при любом h .
Замечание. В общем случае (при произвольной функции F ( x, y ) )
абсолютно устойчивой разностной схеме соответствует разностное уравнение
вида
yi 1 yi
y yi
F ( xi , i 1
) .
h
2
Из такого уравнения не всегда удается выделить y i 1 (для этого может
потребоваться явное решение нелинейного уравнения).
4.5 Задача Коши для обыкновенных дифференциальных уравнений.
Динамическая система
Исследуем поведение пружинного маятника с колеблющееся массой m,
упругостью пружины k и коэффициентом сопротивления (вязкого трения) r в
интервале от 0 до 10 секунд с шагом 0.1 сек при следующих значениях
параметров:
m=1, k=1, r= 0.45
Начальная скорость в положении равновесия равна -0.45
Построим зависимость положения х(t) и скорости x'(t) от времени, а также
фазовую траекторию x (x) и интегральную кривую x(t), x (t).
Используем модификацию метода Эйлера с вычислением производных на
половинном шаге.
Математической моделью данной физической системы является уравнение
x rx kx 0
(4.7)
с начальными условиями
x(0)=0
x (0) 45
Приведем уравнение (4.7) к нормальной системе путем замены
х1=х
х2= x
Получим систему
x1 x2
x 2 kx1 rx2
с начальными условиями
х1(0)=0
х2(0)= -45
Создаем М-файл
function bob
tau=0.1;
t=0:tau:10;
[n1,n2]=size(t);
x=zeros(2,n2);
x(2,1)=-45;
xx=x;
for i=1:n2-1
xx(:,i+1)=x(:,i)+0.5*tau*f(x(:,i));
x(:,i+1)=x(:,i)+tau*f(xx(:,i+1));
end
subplot(2,2,1); plot(t,x(1,:)),xlabel('t'),ylabel('x1'),
grid on
subplot(2,2,2); plot(t,x(2,:)),xlabel('t'),ylabel('x2'),
grid on
subplot(2,2,3);
plot(x(1,:),x(2,:)),xlabel('x1'),ylabel('x2'), grid on
subplot(2,2,4);
plot3(t,x(1,:),x(2,:)),xlabel('t'),ylabel('x1'),
zlabel('x2'), grid on
function u=f(x)
u=[x(2); -x(1)-0.45*x(2)];
end
end
Результат представлен на рисунке 4.6
Рис. 4.6 Задача Коши для мятника
4.6 Средства MATLAB для решения задачи Коши для обыкновенных
дифференциальных уравнений
В
системе
предназначенных
MATLAB
для
имеется
решения
целый
задачи
ряд
Коши
встроенных
для
функций,
обыкновенных
дифференциальных
уравнений.
Краткая
характеристика
этих
функций,
называемых солверами (от английского слова solver, переводимого на русский
язык как "решатель") приведена в таблице 4.1 В них реализованы различные
методы решения задачи Коши для дифференциальных уравнений, поэтому
выбор той или иной функции зависит от характеристик решаемой задачи.
Таблица 4.1
Краткая характеристика солверов
Название функции
ode45
Характеристика функции
реализует одношаговые явные методы Рунге-Кутта 4-го
и 5-го порядка. Это классический метод, рекомендуемый
для начальной пробы решения. Во многих случаях он
дает
хорошие
результаты.
Применяют
часто
эту
функцию, если характеристики задачи неизвестны.
ode23
реализует одношаговые явные методы Рунге-Кутта 2-го
и 4-го порядка. При умеренной жесткости системы ОДУ
и низких требованиях к точности этот метод может дать
выигрыш в скорости решения.
ode113
реализует
многошаговый
метод
Адамса-Башворта-
Мултона переменного порядка. Это адаптивный метод,
который может обеспечить высокую точность решения.
Рекомендуется для решения нежестких задач с высокой
точностью.
ode23tb
реализует неявный метод Рунге-Кутта в начале решения
и
метод,
использующий
формулы
обратного
дифференцирования 2-го порядка в последующем.
Несмотря на относительно небольшую точность, данный
метод может оказаться
например, ode15s.
более эффективным, чем,
ode15s
реализует многошаговый метод переменного порядка (от
1 до 5, по умолчанию 5), использующий формулы
численного
назад".
"дифференцирования
Это
адаптивный метод, его стоит применять, если решатель
ode45
не
обеспечивает
решения
и
система
дифференциальных уравнений жесткая.
ode23s
реализует
одношаговый
метод,
использующий
модифицированную формулу Розенброка 2-го порядка.
Может обеспечить высокую скорость вычислений при
низкой
точности
решения
жесткой
системы
дифференциальных уравнений
ode23t
реализует метод трапеций с интерполяцией. Этот метод
дает
хорошие
описывающих
результаты
при
колебательные
решении
системы
с
задач,
почти
гармоническим выходным сигналом;
bvp4c
служит для решения проблемы граничных значений
систем дифференциальных уравнений вида
y'= f(t,y), F(y(a), y(b), p) = 0 (краевая задача);
pdepe
нужен
для
решения
систем
параболических
и
эллиптических дифференциальных уравнений в частных
производных, введен в ядро системы для поддержки
новых графических функций Open GL, пакет расширения
Partial Differential Equations Toolbox содержит более
мощные средства.
Название этих функций начинается на ode —ordinary differential equetions
(обыкновенные дифференциальные уравнения).
Все решатели, которые представлены в таблице 4.1, могут решать системы
уравнений первого порядка
у' = F(x, y),
причем
для
решения
жестких
систем
уравнений
рекомендуется
использовать лишь специальные функции ode15s, ode23s, ode23t, ode23tb.
Решатели ode15s и ode23t способны найти корни дифференциальноалгебраических уравнений вида
My' = F(x, у),
где М называется матрицей массы.
Все решатели, за исключением
ode23s, который требует постоянства
матрицы массы, могут решать матричного уравнения вида
M(x,y) у' = F(x, у).
Заметим,
что
дифференциальных
ode23tb,
ode23s
уравнений,
дифференциально-алгебраических
ode15s
служат
для
-жестких
уравнений,
решения
жестких
дифференциальных
ode23t
-умеренно
и
жестких
дифференциальных и дифференциально-алгебраических уравнений.
Использование решателей систем ОДУ
В описанных далее функциях для решения систем дифференциальных
уравнений приняты следующие обозначения и правила:
• options — аргумент, создаваемый функцией odeset (еще одна функция —
odeget или bvpget (только для bvp4c)— позволяет вывести параметры,
установленные по умолчанию или с помощью функции odeset или bvpset
соответственно);
• xspan — вектор, определяющий интервал интегрирования [x0 xfinal]. Для
получения решений в конкретных точках разбиения x0, x1,..., xfinal
(расположенные в порядке уменьшения или увеличения) нужно
использовать xspan = [x0 xl ... xfinal];
• у0 — вектор начальных условий;
• p1, р2,… — произвольные параметры, передаваемые в функцию F;
• X, Y — матрица решений Y, где каждая строка соответствует
определенной координате точки разбиения, возвращенном в векторестолбце Х.
Перейдем к описанию функций для решения систем дифференциальных
уравнений:
•
[Х,Y] = solver(@F,хspan,у0) — где вместо solver подставляем имя
конкретного решателя — решается система дифференциальных уравнений
вида у'=F(х,y) на интервале хspan с начальными условиями у0, @F —
дескриптор ODE-функции (также можно задать функцию вида 'F', т.е.
указать имя файл-функции в одинарных кавычках) или, иными словами,
указатель на функцию, которая вычисляет правые части системы
дифференциальных уравнений. Каждая строка в массиве решений Y
соответствует значениям координат точек разбиения, возвращаемым в
векторе-столбце Х;
•
[Х,Y] = solver(@F, хspan, y0. options) — дает решение, подобное
описанному выше, но с параметрами, определяемыми значениями
аргумента options, созданного функцией odeset. Обычно используемые
параметры включают допустимое значение относительной погрешности
RelTol (по умолчанию 1e-3) и вектор допустимых значений абсолютной
погрешности AbsTol (все компоненты по умолчанию равны 1е-6);
•
[Х,Y] = solver(@F,tspan,y0,options,p1,p2,...) — дает решение, подобное
описанному выше, передавая дополнительные параметры p1, р2,... в mфайл F всякий раз, когда он вызывается. Используйте options=[], если
никакие параметры не задаются;
•
[Х,ХE,YE,IE]
=
solver(@F,хspan,y0,options)
—
в
дополнение
к
описанному решению содержит свойства Events, установленные в
структуре options ссылкой на функции событий. Когда эти функции
событий от (х, у) равны нулю, производятся действия в зависимости от
значения трех векторов value, isterminal и direction (их величины можно
установить в m-файлах функций событий). При этом для i-й функции
событий value(i) —значение функции, isterminal(i) — прекратить
интеграцию при достижении функцией нулевого значения; direction(i) = 0,
если все нули функции событий нужно вычислять (по умолчанию),
direction(i) = +1 — только те нули, где функция событий увеличивается,
direction(i) = -1 — только те нули, где функция событий уменьшается.
Выходной аргумент ХЕ — вектор-столбец времен, в которые происходят
события (events), строки YE являются соответствующими решениями, а
индексы в векторе IE определяют, какая из i функций событий (event) равна
нулю при значении независимой переменной, определенный ХЕ. Когда
происходит вызов функции без выходных аргументов, по умолчанию
вызывается выходная функция odeplot для построения вычисленного
решения.
Параметры интегрирования (options) могут быть определены и в m-файле,
и в командной строке с помощью команды odeset. Если параметр определен в
обоих местах, определение в командной строке имеет приоритет.
Решатели используют в списке параметров различные параметры:
• NormControl — управление ошибкой в зависимости от нормы вектора
решения
on
или
off.
Установите
чтобы
'on',
norm(e)<=max(RelTol*norm(y), AbsTol). По умолчанию все решатели
используют более жесткое управление по каждой из составляющих
вектора решения;
• RelTol — относительный порог отбора (положительный скаляр). По
умолчанию равный 1е-3 (т.е. точность оставляет 0.1%) во всех решателях;
оценка ошибки на каждом шаге интеграции e(i) <= max(RelTol*abs(y(i)),
AbsTol(i));
• AbsTol — абсолютная точность (положительный скаляр или вектор, все
компоненты
которого
равны
1е-6).
Скаляр
вводится
для
всех
составляющих вектора решения, а вектор указывает на компоненты
вектора решения. AbsTol по умолчанию принимает значение 1е-6 во всех
решателях;
• Refine - фактор уточнения вывода (положительное целое число) —
умножает число точек вывода на этот множитель. По умолчанию всегда
равен 1, кроме ode45, где он равен четырем. Не применяется, если хspan>2;
• OutputFcn — дескриптор функция вывода (function) — имеет значение в
том случае, если решатель вызывается без явного указания функции
вывода, OutputFcn по умолчанию вызывает функцию odeplot. Эту
установку можно поменять именно здесь;
• OutputSel — индексы отбора (вектор целых чисел) Установите
компоненты, которые поступают в OutputFcn. OutputSel по умолчанию
выводит все компоненты;
• Stats — установка статистики стоимости вычислений. Возможные
значения on или off;
• Jacobian — функция матрицы Якоби (function constant matrix). Данное
свойство следует установить на дескриптор функции FJac (если FJac
возвращает dF/dy) или на имя постоянной матрицы dF/dy;
• Jpattern — график разреженности матрицы Якоби,т.е. имя разреженной
матрицы. Матрица S с S(i,j) = 1, если i-ая составляющая F(х, у) зависит от
j-ой составляющей величины у, и S(i,j) = 0 в противоположном случае;
• Vectorized — векторизованная функция правых частей, возможные
значения которой on или off. Устанавливается на 'on', если функция правых
частей F(x,y) такая, что F(х,[yl y2...]) возвращает вектор [F(х, yl) F(t, y2)
...];
• Events —(function) — ввод дескрипторов функций событий, содержащих
собственно функцию в векторе value, а векторы isterminal и direction (см
выше);
• Mass — матрица массы (constant matrix function), т.е. имя постоянной
матрицы для задач Му' = F(х, у) или дескриптор функции, описывающей
матрицу массы для задач с переменной М вида M(x,y) у' = F(x, у);
• MstateDependence — зависимость матрицы массы от у, возможные
значения: none, weak или strong. Если зависимость отсутствует, то для
уравнений вида Му' = F(х, у) или M(х)*y' = F(х, у) устанавливается 'nоnе'.
'weak' (зависимость слабая) устанавливается для уравнений M(x,y) у' =
F(x, у),
но 'weak' приводит к неявным алгоритмам решения,
использующим аппроксимации при решении алгебраических уравнений.
strong -(зависимость сильная) задается для уравнений M(x,y)у'=F(x,у);
• MassSingular — матрица массы М сингулярная, возможные значения yes
(да), no (нет) или maybe (может быть);
• MvPattern — разреженность (dMv/dy), график разреженности (см
функцию spy) - требуется ввести имя разреженной матрицы S с S(i,j) = 1
для любого k, где (i, k) элемент матрицы массы M(х, у) зависит от
проекции j-ой переменной у, и S(i,j) = 0 в противном случае;
• InitialStep — предлагаемый начальный размер шага, по умолчанию
каждый решатель определяет его автоматически по своему алгоритму;
• InitialSIope — вектор начального уклона ур0, где соответствующий
начальный уклон определяется формулой ур0 = F(х0,y0)/M(х0, y0);
• MaxStep — максимальный шаг, по умолчанию во всех решателях равен
одной десятой интервала хspan;
• BDF (Backward Differentiation Formulas) [on | {off}] — указывает, нужно
ли использовать формулы обратного дифференцирования (методы Gear)
вместо формул численного дифференцирования, используемых в ode15s
по умолчанию. Возможные значения on или off;
• MaxOrder - максимальный порядок ode15s. Возможные значения: [1 | 2 | 3
4 | 5].
Решатели используют в списке различные параметры. В приведенной ниже
таблице они даны для решателей обычных (в том числе жестких)
дифференциальных уравнений.
Таблица 4.2
Параметры решателей, используемые в списке опций.
Параметры
Ode45
Ode23
+
+
+
+
+
+
+
+
+
+
Events
+
+
+
+
+
MaxStep, InitlalStep
+
+
+
+
+
-
-
-
+
+
Mass
-
-
-
+
+
MassConstant
-
-
-
+
-
MaxOrder, BOF
-
-
-
+
-
RelTol,AbsTol
OutputFcaOutputSel,
Refine, Stats
Ode11s
Ode15s
ode23s
Jacobian,
Jpattern, Vectorized
Решатель bvp4c, который мы рассмотрели ранее, имеет очень небольшое
число параметров, но при его использовании можно вводить не только матрицу
Якоби интегрируемой функции, но и матрицу Якоби, содержащую частные
производные функции граничных условий по границам интервала и по
неизвестным параметрам.
Для пояснения понятия жесткости системы дифференциальных уравнений
рассмотрим на ставшем классическом примере — решение уравнения Ван-дерПоля. Уравнение Ван-дер-Поля описывает релаксационные колебания в
электронных устройствах. Это уравнение является жестким. Покажем это.
Уравнение Ван-дер-Поля
u"u k (1 u 2 )u' 0
перепишем в векторной форме:
1
du
u
1 k (1 u12 )
dt
(4.8)
,
u
1
где u u , u1 u, u2 u'
2
2
Величину (1 u1 ) b можно считать порядка единицы, т.е.|b|~1. Находя
собственные значения матрицы (4.8) при k>>1, имеем λ1=bk, λ2=1/(bk).
Число жесткости s определяется как отношение собственных чисел
s=λ1/ λ2=b2k2.
Для k=1000, доходим, что s~106, т.е. система дифференциальных уравнений
является жестким.
В системе MATLAB для нахождения числа жесткости используется
функция cond.
>>A=[0 1;-1 1000];
>>cond(A)
ans=
1e+06
Решим нелинейную систему (4.8) с помощью встроенного в MATLAB
решателя для жестких систем ode23s.
Перед решением нужно записать систему дифференциальных уравнений в
виде ode-функции. Для этого создадим М-файл
function dydt = vdp(t,u)
% формируем вектор столбец
dydt = zeros(2,1);
dydt(1) = u(2);
dydt(2) = 1000*(1 -u(1)^2)*u(2) -u(1);
end
Параметр t будет изменяться в пределах от 0 до 5000, начальные условия
вектора u=[0 1].
Тогда решение решателем ode23s и сопровождающий его график можно
получить, используя следующие команды:
>> [T,Y]=ode23s(@vdp,[0 5000],[0 1]);
>> plot(T,Y(:,1),T,Y(:,2), ...
T(1:(length(T)-1)),diff(T));
>> axis([0,5000,-15,60]);
Решение представлено на рисунке на 4.7
Рис.4.7 Решение уравнения Ван-дер-Поля решателем ode23s
В решателе ode23s происходит автоматическое регулирование шага
интегрирования. Шаг значительно уменьшается вблизи точек с максимальной
производной и, наоборот, вне этих точек существенно увеличивается.
4.7 Решение задачи изгиба консоли функцией MATLAB
Решим рассмотренную выше задачу (4.6) встроенными средствами системы
MATLAB ― функцией ode45
Для рассмотренной задачи текст М-файла в этом случае будет:
function zc
xspan=[0 1];
y0=[0; 0];
ode45(@f,xspan,y0);
function dydx=f(x,y)
c=1.;
dydx=[y(2);M(c*x)*c*(sqrt(1+y(2)^2))^3];
function u=M(x)
u=-1/(sqrt(1+x^2))^3;
end
end
end
Получаем графики прогиба (верхняя кривая на рисунке 4.8) и угла поворота
Рис.4.8 Графики прогиба и угла поворота
Практическая работа 4
Лекция 5
ЗАДАЧА ТЕПЛОПРОВОДНОСТИ
Введение
5.1 Математическая формулировка задачи.
5.2 Метод конечных разностей
5.2.1. Явная схема
5.2.2. Неявная схема
5.2.3. Матричная запись явной схемы
5.2.4. Матричная запись неявной схемы
5.3 Решение задачи теплопроводности на MATLAB
Практическая работа 5
Введение
В задачах расчета сооружений большое значение имеет задача расчета
температуры в конструкции, так как она определяет необходимые условия
функционирования
сооружений,
а
также
термонапряженное
состояние
конструкции.
Пример. Определить распределение температуры по толщине стены для
заданного промежутка времени.
Будем считать, что распределение температуры в стене зависит только от
расстояния точки от наружной части стены и времени (это значит, что стена
имеет достаточно большие размеры). Полагаем, также, что температура на улице
и в помещении в любой момент времени известна. Кроме того, будем считать,
что в исходный момент времени распределение температуры по толщине стены
также известно. Это может быть известно из какого-нибудь предварительного
расчета, например, если исходному моменту предшествовал температурный
период стабильности.
Примечание 1. В действительности по прошествии достаточного времени
распределение температуры мало зависит от начального распределения.
Рис. 5.1
5.1 Математическая формулировка задачи.
u
2u
f ( x, t ) уравнение теплопроводности,
2
t
x
0 x , t 0 пространственно временная область,
(5.1)
u
(
,
t
)
(
t
)
, t 0 краевые условия ,
u( , t ) (t )
u( x,0) ( x ), 0 x начальные условия .
где
x – координата по длине: 0 x ;
t – координата по времени: t 0 ;
u( x, t ) – значение температуры в точке x во время t ;
– коэффициент температуропроводности;
f ( x, t ) – функция источника тепла;
– толщина стены.
Задача (5.1) определена в пространственно-временной области :
{x, t :0 x , t 0}
Рис. 5.2
Отметим, что 0 (0) (0) и ( ) ( ) .
Примечание 2. Поскольку задача (5.1) содержит начальные условия по
времени, то она является задачей Коши.
5.2 Метод конечных разностей
Нанесем на пространственно-временную область прямоугольную сетку.
Пронумеруем узлы сетки:
i 0,1, ... , n – номер узла по длине (ось x );
k 0,1, ... – номер узла по времени ( ось t ).
Таким образом, каждый узел сетки имеет номер, состоящий из двух
компонент (i, k ) . В свою очередь, каждый узел сетки имеет координаты ( xi , t k ) .
Обозначим
uik u ( xi , t k ) ,
– по времени: t k 1 t k ,
h – шаг по длине: h xi 1 xi ,
f i k f ( xi , t k ) ,
0k 0 (t k ), k l (t k ), i ( xi ) .
В зависимости от типа разностной схемы для решения задачи
теплопроводности различаются явная схема и неявная схема, каждая из
которых имеет свои преимущества и недостатки.
5.2.1. Явная схема
Заменим в дифференциальном уравнении задачи (5.1) для внутренних точек
области производные конечно-разностными отношениями. Начальные и краевые
условия будем рассматривать в граничных узлах сетки.
uik 1 uik
i 1,,n 1
uik1 2uik uik1
k
f
,
i
k 0,1,
h2
u k k
0
0
, k 0,1, краевые условия
k
k
u
n
o
ui i , i 0,1,,n начальные условия
(5.2)
В результате конечно-разностной аппроксимации получаем в каждой точке
сетки одно неизвестное u ik и одно условие (граничное или разностное
уравнение).
Первую формулу можно переписать в виде
uik 1 uik
2
(uik1 2uik uik1 ) f i k ,
h
i 1, ..., n 1 .
(5.3)
Таким образом, по этой формуле имея значения uik для всех i , можно
вычислить значения uik 1 для всех i (без решения уравнений).
Рис. 5.3 Явная схема
При решении по явной схеме следует учитывать фактор устойчивости счета,
который накладывает ограничение на величину шага по времени
в
зависимости от величины шага по x :
h2
– условие устойчивости
2
(5.4)
Нарушение условия устойчивости приводит к неправильному результату.
Это ограничение весьма существенно при необходимости решения задачи для
большого периода времени.
Примечание. Условие устойчивости, как правило, получается из
дополнительных аналитических оценок и носит достаточно общий характер для
широкого класса задач. Простой и характерный пример анализа устойчивости,
связанной с выбором шага, рассмотрен в теме “Задача Коши” в данном курсе.
5.2.2. Неявная схема
Существует конечно-разностная аппроксимация, счет по которой не требует
выполнения условия устойчивости.
uik 1 uik
i 1,,n 1
uik11 2uik 1 uik11
k 1
f
,
i
k 0,1,
h2
u k k
0
0
, k 0,1, краевые условия
k
k
u
n
o
ui i , i 0,1,,n начальные условия
(5.5)
Как видно из записи, аппроксимация второй производной по
x
осуществляется для точек (k 1) -го момента времени.
Первую формулу можно переписать в виде
uik 1
h
2
(uik11 2uik 1 uik11 ) uik f i k 1 ,
i 1, ..., n 1
(5.6)
Следовательно, если все величины u ik известны, то величины
u1k 1 , u2k 1 , , unk11
получаем из системы уравнений (5.6) с учетом краевых условий:
u0k 1 0k 1
и
unk 1 lk 1 .
Таким образом, начиная с использования начальных условий
uio i , i 0,1, ..., n ,
последовательно получаем значения неизвестных uik для всех k -ых
моментов времени.
Рис. 5.4 Неявная схема
Примечание. Доказательство устойчивости неявной схемы, как правило,
приводится в математических курсах численных методов. Преимущество
неявной схемы состоит в том, что она устойчива при любом шаге по времени.
Недостатком метода является необходимость решения системы на каждом шаге
по времени.
5.2.3. Матричная запись явной схемы
Явная схема (5.2) может быть представлена в виде
u 0k 1 0k 1 , i 0 ;
uik 1 uik
h
2
(uik1 2uik uik1 ) f i k , i 1, ..., n 1 ;
u nk 1 k 1 , i n .
Введем следующие матрицы и векторы
0
1
E0
,
1
0
u 0k
0k
k
u
1
0
uk , k ,
0
k
u m
k
0 0
1 2 1
A0
,
1
2
1
0 0
f
k
0
0
fk
1
1
, .
f k
n 1
n 1
n
0
Тогда явная схема представима в матричной записи
u k 1 ( E0
h
2
A0 )u k k 1 f k , k 0,1, ... ,
(5.7)
при этом u 0 .
5.2.4. Матричная запись неявной схемы
Неявная схема (5.5) может быть представлена в виде
u 0k 1 0k 1 , i 0 ;
uik 1
h
2
(uik11 2uik 1 uik11 ) uik f i k 1 , i 1, ..., n 1 ;
u nk 1 k 1 , i n .
С учетом принятых обозначений матрично-векторное представление
вышеприведенной системы разностных уравнений будет иметь вид
(E
2
A 0 )u k 1 E0 u k k 1 f k 1
h
здесь E – единичная матрица ( n 1 ) - порядка
(5.8)
1
1
E
1
Обозначим
A E
h2
A0
Тогда матричная запись решения на каждом шаге по времени по неявной
схеме может быть представлена в виде
u k 1 A1 ( E 0 u k k 1 f k 1 )
(5.9)
при этом u 0 .
5.3 Решение задачи теплопроводности на MATLAB
Решим задачу (5.1), т.е. найдем распределение температуры u(x,t) в стене
толщиной
1
и
коэффициентом
температуропроводности
a=0.001.
На
поверхности слева и справа поддерживаются постоянные температуры 110 и 130
соответственно. Начальная температура всюду нулевая.
Применим явную разностную схему, сделав 99 шагов по времени с шагом 5
и разделив стену по толщине на 9 частей. Выведем график распределения
температуры для моментов времени 5, 15, 100, 495.
График распределения температур приведен на рисунке 5.5
Рис.5.5 Распределение температуры для моментов времени 5, 15,100,495
Практическая работа 5
Лекция 6
ЗАДАЧИ ЛИНЕЙНОГО ПРОГРАММИРОВАНИЯ
Введение
6.1 Формулировка задачи линейного программирования
6.2. Пример формулировки и решения задачи линейного
программирования
6.2.1. Геометрическая интерпретация
6.3 Средства MATLAB для решения задач линейного
программирования
Практическая работа 6
Введение
Многие практические задачи, связанные с планированием производства,
экономические задачи, расчет конструкций и т.д. сводятся к, так называемой,
задаче математического программирования. Такая задача состоит в
вычислении максимума (минимума) функции цели Ф(x ) , x ( x1 , x 2 ,..., x n ) при
заданных ограничениях Qi ( x ) 0 ,
i 1,2,...,m ,
m n . Наиболее простым
представителем такой задачи является задача линейного программирования.
6.1 Формулировка задачи линейного программирования
Дана линейная целевая функция
z ( x ) c1 x1 c 2 x 2 ... c n x n
и система m линейных неравенств-ограничений ( m n )
(6.1)
a11 x1 a12 x 2 ... a1n x n b1
a x a x ... a x b
21 1
22 2
2n n
2
,
a m1 x1 a m 2 x 2 ... a mn x n bm
а также дополнительные ограничения
x1 0 , x 2 0 , …, x n 0 .
(6.2)
(6.3)
Требуется найти максимум функции z(x) при выполнении заданных
ограничений (6.2)-(6.3).
Можно записать (6.1) - (6.3) в стандартной и матрично-векторной форме.
Короткая (стандартная) форма
n
z ( x ) ci xi max
(6.1’)
i 1
n
aij x j bi ,
i 1, 2 , ... , m
j 1
x j 0 , j 1, 2 , ... , n
(6.2’)
(6.3’)
Матрично-векторная форма
z(c, x) max
(6.1”)
Ax b
(6.2”)
x0
где
c – вектор целевой функции,
b – вектор ограничений,
A – матрица ограничений.
(6.3”)
6.2. Пример формулировки и решения задачи линейного
программирования
Задачи линейного программирования
достаточно
содержательны с
практической точки зрения и при этом, как правило, имеют решение в классе
точных методов, представителем которых является, например, симплекс-метод.
Рассмотрим практический пример из строительства.
Имеется растворный узел, производящий бетон двух видов (в зависимости
от расхода цемента, песка и щебня). Исходные данные по работе такого узла
могут быть представлены таблицей
Расход сырья на единицу
Вид сырья
продукции
Запас сырья
1й тип бетона
2й тип бетона
Цемент
0.25
0.25
1.5
Песок
0.25
0.5
2.5
Щебень
0.5
0.25
2.5
Цена ед. продукции
2
1
Пусть
x1 – количество выпущенного бетона 1-го типа,
x 2 – количество выпущенного бетона 2-го типа,
тогда
z 2 x1 x 2 – доход растворного узла.
При этом имеет место следующий расход материала:
цемента: 0.25 x1 0.25 x 2 1.5 ,
0.25 x1 0.5 x 2 2.5 ,
песка:
0.5 x1 0.25 x 2 2.5 .
щебня:
(*)
В правой части неравенств (*) присутствует ограничение, связанное с
реальным запасом материалов в растворном узле.
Задача состоит в таком планировании производства (то есть определении
x1 и x 2 – плана выпуска), при котором величина дохода (z ) максимальна, т.е.
найти точку максимума функции z
z 2 x1 x 2
при ограничениях
x1 x 2 6 ,
x 2 x 10 ,
2
1
2 x1 x 2 10 ,
x 0,
1
x 2 0 .
(**)
Здесь неравенства (*) для упрощения расчета умножены на число «4».
Последние два неравенства в (**) являются очевидными, поскольку в нашем
случае объем выпускаемой продукции не может быть отрицательным.
В матричной форме: найти вектор x , при котором функция z (c, x)
достигает максимума, при ограничениях A x b и x 0
где
2
c ,
1
1
A 1
2
1
6
x
2 , b 10 , x 1 .
x2
10
1
6.2.1. Геометрическая интерпретация
Как следует из аналитической геометрии, каждое из неравенствограничений (**) задает некоторую полуплоскость. Их совокупность определяет
некоторый выпуклый многоугольник . Он представлен на рис.6.1. Его
называют многоугольником ограничений. Очевидно, что искомые значения x1
и x 2 принадлежат этому многоугольнику.
Среди точек области необходимо отыскать такую, в которой функция z
принимала бы максимальное значение.
Рис. 6.1.
Выясним геометрический смысл целевой функции z .
Из аналитической геометрии известно, что расстояние d от точки ( x1 , x 2 )
до прямой z 0 , где z ax1 bx 2 c , определяется по формуле
d
a x1 b x 2 c
a 2 b2
1
a 2 b2
z .
Поскольку в рассматриваемой задаче z 0 , то z пропорциональна d и
максимум значения z достигается в точке, максимально удаленной от прямой
z 0.
Теперь задачу можно поставить так.
Среди точек многоугольника найти точку, наиболее удаленную от
прямой
z 2 x1 x 2 0 .
Из рис. 6.1 видно, что такой точкой будет точка пересечения прямых:
x1 x 2 6 ,
2 x1 x 2 10 ,
т.е. одна из вершин многоугольника. Решая эту систему уравнений,
находим координаты искомой точки:
x1 4 , x 2 2 .
Следовательно, в соответствии оптимальному плану (наибольший доход)
нужно выпускать 4 единицы бетона 1-го типа и 2-го единицы бетона 2 типа.
При этом будет получен доход
z 2 4 2 10
Предложенный метод решения – геометрический. Он пригоден для случая
лишь двух-трех переменных.
Перечислим некоторые широко употребляемые термины задач линейного
программирования:
Многогранник ограничений – область допустимых значений.
Допустимое решение – точка, удовлетворяющая всем ограничениям.
Опорное решение – одна из вершин многогранника ограничений. Оно
удовлетворяет всем ограничениям и, по крайней мере, n из них обращает в
равенство.
Оптимальное решение – точка, удовлетворяющая всем ограничениям, в
которой функция цели принимает максимальное значение.
Замечание 1. Из теории задач линейного программирования известно, что
оптимальное решение достигается в одной из вершин многогранника
ограничений. Однако прямой перебор вершин занимает слишком много времени.
Поэтому разработаны методы, которые позволяют вести направленный перебор
вершин (с увеличением величины z). Из таких методов наиболее используемым
является симплекс-метод, алгоритм которого излагается в специальных курсах.
Замечание 2. Многогранник ограничений может быть неограниченным или
сами ограничения могут быть противоречивыми. Эти случаи требуют
специального рассмотрения.
Замечание 3. В реальной задаче линейного программирования помимо
ограничений типа неравенств могут быть ограничения в виде строгих равенств.
Такое равенство можно представить как совокупность двух неравенств.
Например, равенство
a kl x l ... a kn x n bk
может быть записано в виде двух неравенств
a kl x l ... a kn x n bk ,
a kl x l ... a kn x n bk .
6.3 Средства MATLAB для решения задач линейного
программирования
Для решения задач линейного программирования в системе MATLAB
используется функция linprog, соответствующая область поиска для которой
задается следующими условиями:
A*x<=b — линейные неравенства (A-матрица, b- вектор);
Aeq*x=beq - линейные уравнения (Aeq - матрица, beq - вектор);
lb<=x<=ub - ограничения на координаты (lb, ub - векторы);
Целевая функция c'*x в linprog задается вектором коэффициентов c.
Перечислим ниже формы обращения к функции linprog
x= linprog(f,A,b,Aeq,Beq)
x= linprog(f,A,b,Aeq,Beq,lb,ub)
x= linprog(f,A,b,Aeq,Beq,lb,ub,x0)
x= linprog(f,A,b,Aeq,Beq,lb,ub,x0,options)
[x,fval]= linprog(. . .)
[x,fval,exitflag]= linprog(. . .)
[x,fval,exitflag,output]= linprog(. . .)
[x,fval,exitflag,output,lambda]= linprog(. . .)
где
х — искомый вектор;
с — вектор коэффициентов функции цели;
fval — искомое минимальное значение функции цели;
exitflag — параметр, характеризующий вычислительный процесс,( если его
значение больше нуля, то результат найден с требуемой точностью, нуль достигнуто максимальное значение итераций, меньше нуля - решение не
найдено);
output — структура, содержащая сведения о процессе оптимизации;
lambda — структура, содержащая сведения о множителях Лагранжа для
решения х.
х0— стартовая (начальная) точка.
Отсутствующие ограничения в списке параметров задаются квадратными
скобками.
В функции linprog по умолчанию используется так называемый прямодвойственный алгоритм, при использовании которого одновременно решается
как исходная задача линейного программирования, так и двойственная к ней
(при этом решение двойственной задачи не выдается).
Решим задачу линейного программирования (**) с использованием
функции linprog
Текст М-файла:
Результаты расчета
Введите n= 2
Введите m= 3
Вектор коэффициентов функции цели, с
2.0000
1.0000
Матрица коэффициентов ограничений, A
1.0000
1.0000
1.0000
2.0000
2.0000
1.0000
Вектор правых частей ограничений, b
6.0000
10.0000
10.0000
Optimization terminated.
Вектор решения, x
4.0000
2.0000
Практическая работа 6
