Учебное пособие 646

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

>> A=[1 2 3 4 5; 1 3 5 7 9; 7 6 8 1 2] A =

>>B=[1; 2; 3] B =

1

2

3

>>X=A\B

X = -0.8571

0

1.2143

0

-0.3571

Сделаем проверку:

>> A*X ans =

1.0000

2.0000

3.0000

>> A=[7 -5 -2 -4; -3 2 1 2; 2 -1 -1 -2; -1 0 1 2; 0 -1 1 2] A =

>>B=[8; -3; 1; 1; 3] B =

8 -3 1 1 3

>>X=A\B

Warning: Rank deficient, rank = 2 tol = 8.8121e-015.

X=

2.0000

0

0

1.5000

Проверка:

>> A*X ans =

8.0000 -3.0000 1.0000

1.0000

3.0000

Задание № 2 к выполнению лабораторной работы

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

Контрольные вопросы

1.В чем состоит метод решения системы линейных уравнений с помощью обратной матрицы?

2.В чем состоит метод Гаусса решения системы линейных уравнений?

3.В чем заключаются преимущества и недостатки этих методов?

4.Какие системы называются недоопределенными и переопределенными?

5.Каким методом решаются такие системы?

6.Что такое базисное решение?

7.Как MATLAB решает несовместные системы?

ЛАБОРАТОРНАЯ РАБОТА № 4

ЧИСЛЕННОЕ ДИФФЕРЕНЦИРОВАНИЕ

Численное дифференцирование реализуется следующей функцией:

- diff X – возвращает конечные разности смежных элементов массива X . Если X – вектор, то diff X возвращает вектор разностей соседних элементов X i 1 X i , у которого количество элементов на единицу меньше чем у исходного вектора X . Если X – матрица, то diff X возвращает матрицу разностей соседних строк;

- diff X , n – возвращает конечные разности порядка n смежных элементов массива X ;

- diff X , n, dim – возвращает конечные разности порядка n .

Если dim 1, то определяется разность соседних строк, если dim 2 , то определяется разность соседних столбцов, если dim 3 , то определяется разность соседних элементов в третьей размерности матрицы X . Если величина dim превышает размерность матрицы, то выдается сообщение об ошибке.

Примеры:

по x . Как и следовало ожидать, полученный график (рис. 4) очень близок к графику функции cos x sin x .

Задание к выполнению лабораторной работы

На собственных примерах разобрать работу функции diff с векторами и матрицами. Построить график производ-

ной функции f x взятой из своего варианта к лабораторной работе № 3 на отрезке a, b .

Контрольные вопросы

1.Запишите простейшие формулы для численного вычисления производных первого и второго порядка.

2.Как получить эти формулы из формулы Тейлора?

3.Каков порядок точности этих формул?

4.Дайте геометрическую интерпретацию формул для производной первого порядка.

5.С помощью какой функции производится численное дифференцирование в системе MATLAB?

6.Какие величины вычисляет эта функция и ее разновидности, если X – вектор?

7.Какие величины вычисляет эта функция и ее разновидности, если X – матрица?

8.Как производится построение графика производной?

9.Почему в примере построения графика производной делалось повторное определение массива X ?

ЛАБОРАТОРНАЯ РАБОТА № 5

ЧИСЛЕННОЕ ИНТЕГРИРОВАНИЕ

Численное интегрирование заключается в приближенном вычислении определенного интеграла вида

24

b

f x dx

a

одним из многочисленных численных методов. Приведенные ниже функции выполняют численное ин-

тегрирование методом трапеций и методом трапеций с накоп-

лением:

- trapz Y – возвращает определенный интеграл, используя

интегрирование методом трапеций с единичным шагом между отсчетами. Если Y – вектор, то функция возвращает интеграл элементов вектора Y . Если Y – матрица, то функция возвращает вектор-строку, содержащую интегралы каждого столбца этой матрицы;

- trapz X ,Y – возвращает определенный интеграл от функ-

ции Y по переменной X , используя метод трапеций (пределы интегрирования в этом случае задаются начальным и конечным элементами вектора X ).

Примеры:

>>Y=[1 2 3 4]

>>trapz(Y) ans =

7.5000

>>X=0:pi/50:pi/2;

>>Y=cos(X);

>>trapz(X,Y)

ans = 0.9997

- cumtrapz Y – возвращает численное значение определенно-

го интеграла для функции, заданной ординатами в векторе или матрице Y с шагом интегрирования, равным единице (интегрирование методом трапеций с накоплением), В случае когда шаг отличен от единицы, но постоянен, вычисленный интеграл достаточно умножить на величину шага. Для векторов эта

25

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

- cumtrapz X ,Y – выполняет интегрирование с накоплением

от функции Y по переменной X , используя метод трапеций (пределы интегрирования в этом случае задаются начальным и конечным элементами вектора X ).

Пример:

>>Y=[1 2 3 4];

>>cumtrapz(Y)

ans =

0 1.5000 4.0000 7.5000

Приведенные ниже функции осуществляют интегрирование, используя квадратурную формулу Симпсона или метод Гаусса-Лобатто. Функция quad выполняет интегрирование по

методу низкого порядка, используя рекурсивное правило Симпсона. Но она может быть более эффективной при негладких подынтегральных функциях или при низкой требуемой точности вычислений. Алгоритм quad в MATLAB-6 изменен

по сравнению с предшествующими версиями, точность по умолчанию по сравнению с версиями 5.* повышена в 1000 раз

и составляет 10 6 Новая функция quadl (квадратура Лобатто)

использует адаптивное правило квадратуры Гаусса-Лобатто очень высокого порядка. Поэтому достижимая точность интегрирования гладких функций в MATLAB-6 также значительно выше, чем в предшествующих версиях

-quad @ fun, a,b – возвращает численное значение определенного интеграла от заданной функции на отрезке a, b ;

-quad @ fun, a,b,tol – возвращает численное значение определенного интеграла от заданной функции на отрезке a, b с

26

заданной относительной погрешностью tol . Можно также использовать вектор tol rel _ tol abs _ tol , чтобы точно опре-

делить комбинацию относительной и абсолютной погрешности.

Примеры:

>>quad('sin(x)',0,pi) ans =

2.0000

>>quadl('sin(x)',0,pi) ans =

2.0000

Задание к выполнению лабораторной работы

Вычислить интеграл из своего варианта с использованием всех приведенных выше функций.

Варианты исходных данных

Контрольные вопросы

1.В чем состоит метод трапеций вычисления определенного интеграла?

2.В чем состоит метод Симпсона вычисления определенного интеграла?

3.В чем состоит метод Гаусса-Лобатто вычисления определенного интеграла?

28

4.Какие еще методы вычисления определенного интеграла вы знаете и в чем они заключаются?

5.Какая функция осуществляет метод трапеций?

6.Какая функция осуществляет метод Симпсона?

7.Какая функция осуществляет метод Гаусса-Лобатто?

8.Какой метод вычисления интегралов из рассмотренных выше точнее и почему?

ЛАБОРАТОРНАЯ РАБОТА № 6

ЗАДАЧА КОШИ ДЛЯ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ

Анализ поведения многих систем и устройств в динамике, а также математическое описание природных процессов обычно базируется на решении обыкновенных дифференци-

альных уравнений (ОДУ) и их систем. Их, как правило, пред-

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

dydx y f x, y

с граничными условиями y x0 y0 . Решение ищется на от-

резке x0 , b .

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

-ode45 – одношаговые явные методы Рунге-Кутта 4-го и 5-го порядка. Это классический метод, рекомендуемый для начальной пробы решения. Во многих случаях он дает хорошие результаты;

-ode23 – одношаговые явные методы Рунге-Кутта 2-го и 3-го порядка. При умеренной жесткости системы ОДУ и низких

29

требованиях к точности этот метод может дать выигрыш в скорости решения;

-ode113 – многошаговый метод Адамса-Башфорта-Мултона переменного порядка. Это адаптивный метод, который может обеспечить высокую точность решения. Решатель ode113 оказывается особенно эффективным для нежестких систем дифференциальных уравнений, правые части которых вычисляются по сложным формулам.

Если все попытки применения рассмотренных выше решателей не приводят к успеху, то, возможно, что решаемая система является жесткой. Для решения жестких систем рекомендуется использовать только специальные решатели:

-ode15s – многошаговый метод Гира переменного порядка (от 1 до 5, по умолчанию 5). Это адаптивный метод, его следует применять, если решатель ode45 не обеспечивает решения;

-ode23s – одношаговый метод, использующий модифицированную формулу Розенброка 2-го порядка. Может обеспечить высокую скорость вычислений при низкой точности решения жесткой системы дифференциальных уравнений;

-ode23t – метод трапеций с интерполяцией. Этот метод дает хорошие результаты при решении задач, описывающих колебательные системы с почти гармоническим выходным сигналом;

-ode23tb – неявный метод Рунге-Кутта в начале решения и метод, использующий формулы обратного дифференцирования 2-го порядка в последующем. Несмотря на сравнительно низкую точность, этот метод может оказаться более эффективным, чем ode15s .

Спомощью команды odeset можно определить следую-

щие параметры интегрирования (приводим не все):

-RelTol – относительная точность решения (скаляр или вектор) Скаляр вводится для всех составляющих вектора решения, а вектор для каждой компоненты вектора решения. По

30

умолчанию во всех решателях относительная точность реше-

ния равна 10 3 ;

- AbsTol – абсолютная точность решения (скаляр или вектор) Скаляр вводится для всех составляющих вектора решения, а вектор для каждой компоненты вектора решения. По умолча-

нию во всех решателях абсолютная точность равна 10 6 ;

- Refine – фактор уточнения вывода – умножает число точек

вывода на этот множитель. По умолчанию всегда равен 1, кроме ode45 , где он равен 4. Не применяется, если интервал интегрирования больше 2;

- InitialStep – предлагаемый начальный размер шага, по умол-

чанию каждый решатель определяет его автоматически по своему алгоритму;

- MaxStep – максимальный шаг, по умолчанию во всех реша-

телях равен одной десятой интервала интегрирования;

- MaxOrder – максимальный порядок, используемый решате-

лем ode15s (1, 2, 3, 4, 5).

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

>> options=odeset('RelTol',[1e-7 1e-8],'AbsTol',3e-5);

Функция odeset без параметров возвращает все имена параметров и их допустимые значения:

>> odeset

AbsTol: [ positive scalar or vector {1e-6} ] RelTol: [ positive scalar {1e-3} ]

NormControl: [ on | {off} ] OutputFcn: [ function ] OutputSel: [ vector of integers ]

Refine: [ positive integer ] Stats: [ on | {off} ]

InitialStep: [ positive scalar ] MaxStep: [ positive scalar ] BDF: [ on | {off} ] MaxOrder: [ 1 | 2 | 3 | 4 | {5} ]

Jacobian: [ matrix | function ] JPattern: [ sparse matrix ]

Vectorized: [ on | {off} ] Mass: [ matrix | function ]

MtateDependence: [ none | weak | strong ] MvPattern: [ sparse matrix ]

MassSingular: [ yes | no | {maybe} ] InitialSlope: [ vector ]

Events: [ function ]

Перейдем к описанию функций для решения систем дифференциальных уравнений:

- X ,Y solver @ f ,tspan, y0 – интегрирует систему ОДУ на отрезке tspan . Здесь Y - матрица решений, где каждая

строка соответствует значению аргумента, возвращаемому в векторе-столбце X ; слово solver означает имя одного из решателей; @ f - дескриптор правых частей системы; y0 - век-

тор начальных условий;

- X ,Y solver @ f ,tspan, y0,options – дает решение, по-

добное описанному выше, но с параметрами интегрирования, определяемыми значениями аргумента options , созданного

функцией odeset ;

- X ,Y solver @ f ,tspan, y0, options, p1, p2, – дает реше-

ние, подобное описанному выше, передавая дополнительные параметры p1 , p2 , в m -файл f всякий раз, когда он вы-

зывается. Если никакие параметры интегрирования не задаются, то вместо options следует использовать [] .

Рис. 5 Приведем пример решения дифференциального уравне-

ния y x cos y с начальным условием y 1,8 2,6 на отрезке 1,8; 2,8 .Перед решением необходимо записать функцию в

правой части уравнения в виде m -файла. Для этого в главном меню выберем File►New►M-file и введем

function f=fun(x,y) f=x+cos(y)

Сохраним этот m -файл и закроем окно. Затем получаем решение y x и строим его график (рис. 5):

>>[X,Y]=ode45('fun',[1.8 2.8],2.6);

>>plot(X,Y)

Приведем пример решения системы ОДУ. В качестве объекта исследования возьмем модель Лотка-Вольтерра борь-

бы за существование. Обозначим: y1 t - число жертв, y2 t -

число хищников. Число хищников и жертв в течение времени t изменяется по закону

y1 Py1 py1y2

y2 Ry2 ry1y2

где Py1 - увеличение числа жертв в отсутствие хищников,Ry2 - уменьшение числа хищников в отсутствие жертв. Вероятность поедания хищником жертвы пропорциональна числу y1y2 , при этом слагаемое py1y2 соответствует вымиранию жертв, а ry1y2 - появлению хищников. Известно, что решением этой системы на плоскости y1y2 является замкнутая

кривая. Возьмем для примера значения коэффициентов P 3 , R 2 , p r 1, и считаем, что в начальный момент времени

было три жертвы и четыре хищника. Таким образом, имеем математическую задачу – решить систему ОДУ

y1 3y1 y1y2y2 2y2 y1y2

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

y1 0 3 , y2 0 4 .

Решим эту систему ОДУ на отрезке t 0, 100 при по-

мощи решателя ode45 :

Создаем m -файл правых частей ОДУ function f=LotVol(t,y)

f=[3*y(1)-y(1)*y(2); -2*y(2)+y(1)*y(2)]

Задаем вектор начальных условий

>> y0=[3; 4];

Вызов решателя

>> [T,Y]=ode45('LotVol',[0 100],y0);

Вывод графика решения системы ОДУ в виде параметрической кривой

>> subplot(1,1,1); plot(Y(:,1),Y(:,2));

Вывод пояснений на график

>> xlabel('жертвы'); ylabel('хищники')

35

Графики приближенных решений, приведенные на рис. 6 и рис. 7 отличаются друг от друга. Приближенное решение, полученное решателем ode23s , намного точнее, чем решателем ode45 . Дело в том, что уравнения Лотка-Вольтерра являются примером так называемых жестких уравнений, для решения которых следует использовать специально приспособленные решатели.

Задание № 1 к выполнению лабораторных работ

Найти решение уравнения y f x, y с начальным условием y x0 y0 на отрезке x0 , b и построить его график.

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

36

Варианты исходных данных

при начальных условиях

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

37

Варианты исходных данных