Вычисление последующего элемента через предыдущие в матлаб. Простые вычисления в MATLAB


1. Command Window (Окно команд).

Математические выражения пишутся в командной строке после знака приглашения « >> ». Например,

Для выполнения действия нажмем клавишу «Enter».

По умолчанию программа записывает результат в специальную переменную ans.

Для сохранения результата под другим именем используют имя переменной и знак присваивания « = », например

>> z = 1.25 /3.11

Редактировать в Command Window можно только текущую командную строку. Для того чтобы отредактировать ранее введенную команду, необходимо установить курсор в строку ввода и воспользоваться клавишами « » или« ».

Если команда заканчивается «;», то результат её действия не отображается в командной строке.

Командное окно можно закрыть кнопкой « », а кнопка « » служит для отделения окна от интерфейса системы. Вернуться к стандартной форме окна можно с помощью меню:

Главное Меню Desktop Desktop Layout Default .

Очистить командное окно можно с помощью меню:

Главное Меню Edit Clear Command Window .

Главное меню системы MatLab.

Главное меню MatLab содержит следующие пункты:

File (Файл) – работа с файлами;

Edit (Правка) – редактирование;

View (Вид) – управление окнами;

Web – связь с фирмой – разработчиком через Интернет;

Window (Окно) – связь с окнами системы;

Help (Справка) – связь со справочной системойMatLab.

Панель инструментов системы MATLAB.

Кнопки панели инструментов имеют следующие назначения:

New file (Создать) – выводит окна редакторов файлов;

Open file (Открыть) – открывает окна загрузки файлов;

Cut (Вырезать) – вырезает выделенный фрагмент и помещает в буфер обмена;

Copy (Копировать) – копирует выделенный фрагмент в буфера обмена;

Paste (Вставить) – переносит выделенный фрагмент из буфера обмена в строку ввода;

Undo (Отменить) – отменяет результата предыдущей операции;

Redo (Повторить) – повторяет результат предыдущей операции;

Simulink – создает модель Simulink (расширения MatLab );

Help Window (Помощь) – открывает окна справки.

4. Формат вывода результата вычислений.



При вводе вещественных чисел для отделения дробной части используется точка !

>> s = 0.3467819

Результат вычислений выводится в формате short (краткая запись числа), который определяется по умолчанию. Можно поменять формат на long (длинная запись числа).

>> format long

0.34678190000000

В списке Numerical Format имеются форматы чисел

short – краткая запись числа;

long – длинная запись числа;

short e – краткая запись числа в формате с плавающей точкой;

long e – длинная запись числа в формате с плавающей точкой;

short g – вторая форма краткой записи числа;

long g – вторая форма длинной записи числа;

Формат отображения числовых данных можно установить в меню File (файл) пункт Preferences (предпочтения). Перейти на вкладку Command Window (окно команд). В опции Text display (отображение текста), в списке Numeric format (числовой формат) установить short g , в опции Numeric display (отображение чисел) установить compact . Эти форматы вывода соответствуют выводу чисел в универсальной форме из пяти значащих цифр и с подавлением пробела между строками.

Основы вычислений в MatLab.

Для выполнения простейших арифметических операций в MatLab применяются операторы:

· сложение и вычитание +, – ;

· умножение и деление *, / ;

· возведение в степень ^ .

Некоторые специальные переменные :

ans – результат последней операции без знака присваивания;

eps – относительная погрешность при вычислениях с плавающей точкой;

pi – число ;

i или j – мнимая единица;

Inf – бесконечность ;

NaN – неопределенное значение.

Некоторые встроенные элементарные функции MatLab :

exp(x) – экспонента числа x;

log(x) – натуральный логарифм числа x;

sqrt(x) – квадратный корень из числа x;

abs(x) – модуль числа x;

sin(x), cos(x), tan(x), cot(x) – синус, косинус, тангенс, котангенс числа x;

asin(x), acos(x), atan(x), acot(x) – арксинус, арккосинус, арктангенс, арккотангенс числа x;

sec(x), csc(x) – секанс, косеканс числа x;

round(x) – округление числа x до ближайшего целого;

mod(x,y) – остаток от целочисленного деления x на y с учетом знака;

sign(x) – возвращение знака числа x.

Вычислим значение выражение

>> exp(–2.5)*log(11.3)^0.3 – sqrt((sin(2.45*pi)+cos(3.78*pi))/tan(3.3))

Если оператор не удается разместить в одной строке, то возможно продолжение его ввода в следующей строке, если в конце первой строки указать знак продолжения «…», например,

>> exp(–2.5)*log(11.3)^0.3 – ...

sqrt((sin(2.45*pi)+cos(3.78*pi))/tan(3.3))

Функции для работы с комплексными числами :

Ввод комплексного числа

>> z = 3 + 4i

3.0000 + 4.0000i

Функции abs(z), angle(z) возвращают модуль и аргумент комплексного числа , где , ;

complex(a,b) конструирует комплексное число по его действительной и мнимой части:

conj(z)возвращает комплексно-сопряженное число;

imag(z), real(z) возвращает мнимую и действительную часть комплексного числа z.

6. Векторы .

Ввод, сложение, вычитание, умножение на число.

Вектор в MatLab формируется с помощью оператора квадратные скобки . При этом элементы вектора-столбца разделяют точкой с запятой «;», а элементы вектора-строки разделяют пробелом « » или запятой « , ».

Введем вектор-столбец .

>> x =

Введем вектор-строку .

>> y =

Для транспонирования вектора применяют апостроф «’ »:

>> z = y’

Для нахождения суммы и разности векторов используются знаки « + » и «– »:

>> с = x + z

Умножение вектора на число осуществляется как справа, так и слева при помощи знака « * ».

Векторы могут быть аргументами встроенных функций, например,

>> d = sin(c)

Для обращения к элементам векторов используются скобки (), например,

>> x_2 = x(2)

Последний элемент вектора можно выбрать, набрав команду

>> X_end = x(end)

Из нескольких векторов можно составить один, например

>> r =

1.3 5.4 6.9 7.1 3.5 8.2

Символ двоеточие « : » используется для выделения нескольких элементов из вектора, например

>> w = r(3:5)

Символ двоеточие « : » также позволяет заменять элементы вектора, например,

>> r(3:5)= 0

1.3 5.4 0 0 0 8.2

Символ « : » также можно использовать для построения вектора, каждый элемент которого отличается от предыдущего на постоянное число, т.е. шаг, например

>> h =

1 1.2 1.4 1.6 1.8 2

Шаг может быть отрицательным (в этом случае начальное число должно быть больше конечного).

Шаг, равный единице, можно не указывать

>> k =

Основные функции для работы с векторами .

  • length(x) – определение длины вектора x;
  • prod(x)– перемножение всех элементов вектора x;
  • sum(x) – суммирование всех элементов вектора x;
  • max(x)– нахождение максимального элемента вектора x;
  • min(x)– нахождение минимального элемента вектора x.

Если вызвать функцию min или max с двумя выходными аргументами = min(x),

то первой переменной присваивается значение минимального (максимального элемента), а второй переменной присваивается номер этого элемента.

7 Матрицы .

Различные способы ввода матрицы .

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

>> A =

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

>> A = .

Алгоритм, реализованный в операции \ состоит из следующих шагов:

  • 1. В тривиальном случае, если A разреженная и диагональная (разреженным матрицам в MATLAB посвящен раздел Разреженные матрицы), решение находится по простой формуле x k = b k /a kk , где k=1,2,...n.
  • 2. Если A - квадратная, разреженная и ленточная матрица, то используется солвер для ленточных матриц. При этом могут быть два варианта:

    a. Если A - трехдиагональная матрица и b - один столбец из вещественных чисел, тогда система решается гауссовым исключением (его операции совершаются только с элементами на диагоналях). Если в процессе исключения для сохранения устойчивости потребуется сделать перестановки строк, или если матрица A не является трехдиагональной, то работает следующий пункт.
    b. Если A - ленточная с плотностью ненулевых элементов больше параметра bandden , по умолчанию равного 0.5, или не выполняются условия предыдущего пункта, тогда, в зависимости от типа A и b (double или single ) вызываются следующие процедуры библиотеки LAPACK:

    A и b вещественные типа double - процедуры DGBTRF, DGBTRS
    A и b комплексные типа double - процедуры ZGBTRF, ZGBTRS
    A или b вещественные типа single - процедуры SGBTRF, SGBTRS
    A или b комплексные типа single - процедуры CGBTRF, CGBTRS

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

    nz = 1 (на диаг.№-2) + 6 (на диаг.№-1) + 7 (на диаг.№0) + 6 (на диаг.№1) + 1 (на диаг.№2) + 1 (на диаг.№3) = 22,

    а число всех элементов в ленте

    band = 5 (на диаг.№-2) + 6 (на диаг.№-1) + 7 (на диаг.№0) + 6 (на диаг.№1) + 5 (на диаг.№2) + 4 (на диаг.№3) = 33

    и плотность ненулевых элементов будет 2/3. Поскольку 2/3 > 0.5, то будет применен солвер для ленточных матриц
    Параметр bandden , как и другие параметры, управляющие алгоритмами для разреженных матриц в MATLAB, устанавливается при помощи функции spparms.

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

  • 3. Если A - верхняя или нижняя треугольная матрица, то применяется метод обратной подстановки, или, соответственно, метод прямой подстановки, в котором из последнего (или первого уравнения) находится компонента решения и далее найденные компоненты подставляются в следующие уравнения для нахождения очередных компонент решения системы уравнений.
  • 4. Если A - может быть приведена перестановками к треугольной форме, то выполняется соответствующее приведение, а далее система уравнений решается как и в п.3.
  • 5. Если A - симметричная матрица (или, в комплексном случае, эрмитова), на ее диагонали находятся только положительные вещественные элементы, то в зависимости от того, разрежена A или нет, выполняется п. а) или п. б).

    a. Если А не является разреженной, тогда делается попытка выполнить разложение Холецкого A = R T R с последующим решением систем с треугольными матрицами R T и R: R T y = b и Rx = y. При этом вызываются следующие процедуры библиотеки LAPACK:

    A и b вещественные и типа double - DLANGE, DPOTRF, DPOTRS, DPOCON
    A и b комплексные и типа double - ZLANGE, ZPOTRF, ZPOTRS, ZPOCON
    A и b вещественные и типа single - SLANGE, SPOTRF, SPOTRS, SDPOCON
    A и b комплексные и типа single - CLANGE, CPOTRF, CPOTRS, CPOCON

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

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

    p = symmmd(A) - в вектор p записаны перестановки

    R = chol(A(p, p));

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

    и вторая, с множителем разложения и с занесением компонент решения на соответствующие позиции в векторе x
    x(p, :) = R \ y

  • 6. Если A - хессенбергова матрица, то она приводится к верхней треугольной (с соответствующей модификацией правой части) и затем система решается подстановками
  • 7. Если A - квадратная матрица, не удовлетворяющая условиям пунктов 1-6, то в зависимости от того, разрежена она или нет, выполняются следующие действия

    a. Если A не является разреженной матрицей, то при помощи исключения Гаусса с выбором ведущего элемента (для обеспечения устойчивости разложения) выполняется LU-разложение матрицы A = LU, где

    U - верхняя треугольная матрица
    L - матрица, перестановками строк сводящаяся к треугольной

    и решение системы Ax = b находится последовательным решением систем с треугольными матрицами Ly = b, Ux = y.
    Для выполнения LU-разложения вызываются следующие процедуры библиотеки LAPACK:

    A и b вещественные и типа double - DLANGE, DGESV, DGECON
    A и b комплексные и типа double - ZLANGE, ZGESV, ZGECON
    A и b вещественные и типа single - SLANGE, SGESV, SGECON
    A и b комплексные и типа single - CLANGE, CGESV, CGECON

    b. Если A - разреженная матрица, то производится перестановка столбцов с целью уменьшения заполнения множителей L и U в процессе нахождения разложения. Далее перестановкой строк в процессе LU-разложения обеспечивается вычислительная устойчивость, после которого снова решаются системы с треугольными матрицами. Для выполнения LU-разложения разреженной матрицы используются процедуры библиотеки UMFPACK

    8. Если предыдущие пункты алгоритма не сработали и, следовательно, A не является квадратной матрицей, то все определяется ее разреженностью и работает один из приведенных ниже пунктов:

    a. Если A не является разреженной матрицей, то выполняется QR-разложение AP = QR, где P - матрица перестановки столбцов, Q - ортогональная матрица (Q T Q = I), а R - верхняя треугольная. Если A размера m на n, то Q размера m на m, а R размера m на n. Далее решение системы находится так:

    x = P*(R \ (Q" * b))

    поскольку из Ax = b и AP = QR следует: QRx = bP и Rx = Q T bP.

    b. В случае разреженной и прямоугольной матрицы A нет смысла вычислять матрицу Q, поскольку она скорее всего окажется заполненной. Поэтому в ходе алгоритма QR-разложения вычисляется c = Q T b (т.е. модифицированная правая часть) и решается система Rx = c с треугольной матрицей для получения решения исходной системы Ax = b.

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

Функция linsolve для решения систем линейных уравнений

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

Функция linsolve, вызванная в самом простом варианте с двумя входными аргументами и одним выходным аргументом
x = linsolve(A, b)
решает систему Ax = b одним из способов, в зависимости от того, квадратная матрица, или нет.

  • Если A - квадратная матрица, то предварительно вычисляется ее LU-разложение и затем решаются две системы с треугольными матрицами L и U.
  • Если A - прямоугольная матрица, то предварительно вычисляется ее QR-разложение и затем решается система с треугольными матрицами.

(Более подробно про действия с полученными множителями разложения написано в разделе ). При этом, если матрица A плохообусловлена или A - матрица неполного ранга, то выводится соответствующее предупреждение, например:

A = ; b = ; x = linsolve(A,b) Warning: Rank deficient, rank = 1, tol = 4.4686e-015. x = 0 0 2.0000

Для подавления вывода данного сообщения в командное окно следует вызывать функцию linsolve с двумя выходными аргументами
= linsolve(A, b)
тогда в x запишется полученное решение, а в r - или ранг матрицы (если A является прямоугольной матрицей), или обратная величина к оценке для ее числа обусловленности (см. раздел ), например

A = ; b = ; = linsolve(A,b) x = -1.0000e+000 9.9999e-001 3.3307e-006 r = 6.9444e-013

Основные преимущества функции linsolve над операцией \ проявляются при указании типа матрицы системы уравнений. Для этого перед вызовом функции linsolve надо заполнить структуру с информацией о матрице со следующими полями

  • SYM - симметричная;
  • LT - нижняя треугольная;
  • UT - верхняя треугольная;
  • UHESS - хессенбергова;
  • POSDEF - симметричная, положительно определенная;
  • RECT - прямоугольная;
  • TRANSA - надо ли решать систему с матрицей, транспонированной к заданной.

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

LT UT UHESS SYM POSDEF RECT TRANSA
true false false false false true/false true/false
false true false false false true/false true/false
false false true false false false true/false
false false false true true false true/false
false false false false false true/false true/false

Если матрица системы положительно определена, то это обязательно учесть при решении, поскольку для положительно определенных матриц решение основано на разложении Холецкого, которое требует меньше операций, чем LU-разложение, используемое при решении систем с квадратными матрицами общего вида. В этом несложно убедиться при помощи следующего примера, в котором создается симметричная положительно определенная матрица (матрица из случайных чисел складывается с транспонированной к ней и на диагональ добавляются большие числа) и система уравнений с этой матрицей сначала решается как система с матрицей общего вида (opts.SYM = false и opts.POSDEF = false), а затем - как с симметричной и положительно определенной матрицей (opts.SYM = true и opts.POSDEF = true).

% задаем все поля структуры ops, кроме SYM и POSDEF opts.TRANSA = false; opts.UHESS = false; opts.RECT = false; opts.UT = false; opts.LT = false; % создаем вектор размеров матриц N = 2.^(8:12); % создаем пустые массивы для записи времен решения TimeGeneral = ; TimePosSym = ; % в цикле создаем матрицы и сравниваем времена решения for n = N % создаем симметричную и положительно определенную матрицу % и вектор правой части A = rand(n); A = A + A" + 2*n*eye(n); b = sum(A, 2); % решаем систему как систему с матрицей общего вида opts.SYM = false; opts.POSDEF = false; Tstart = cputime; x = linsolve(A,b, opts); Tend = cputime; TimeGeneral = ; % решаем систему как систему с симметр и полож. опред матрицей opts.SYM = true; opts.POSDEF = true; Tstart = cputime; x = linsolve(A,b, opts); Tend = cputime; TimePosSym = ; end % выводим графики временных затрат figure loglog(N, TimeGeneral, N, TimePosSym) legend("TimeGeneral", "TimePosSym")

Получающиеся вычислительные затраты этих способов решения системы приведены на графике ниже

Разумеется, если матрица треугольная, то очень важно это указать, поскольку решение системы с треугольной матрицей выполняется за O(n 2) операций, а если к системе с треугольной матрицей применить алгоритм решения, предназначенный для матрицы общего вида, то понадобится порядка O(n 3) операций.

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

A = ; b = ; задать true для поля UT (а все остальные установить в false) opts.UT = true; opts.TRANSA = false; opts.LT = false; opts.UHESS = false; opts.SYM = false; opts.POSDEF = false; opts.RECT = false; то при решении системы функция linsolve будет рассматривать матрицу как верхнюю треугольную и выберет нужные ей элементы из верхнего треугольника A x = linsolve(A,b, opts) При этом получится решение x = 1 1 1 соответствующее матрице A = ;

Влияние обусловленности матрицы на точность решения системы с ней

В этом разделе мы рассмотрим, как ошибки в элементах вектора правой части и в матрице системы линейных уравнений Ax = b могут повлиять на решение этой системы. Обратимся сначала к случаю внесения возмущений в вектор правой части b. Итак, мы решаем две системы

причем предполагается, что каждую из систем мы решаем точно, а основной возникающий вопрос - насколько будет отличаться решение x системы (1) от решения системы (2) с возмущением в правой части δb. Это довольно важный вопрос, поскольку элементы правой части могут быть получены как результат некоторых измерений, соответственно содержащий погрешность δb k в каждой компоненте b k . Хотелось бы, чтобы при измерении одной и той же величины (каждый раз со своими, небольшими погрешностями) соответствующие решения систем с мало отличающимися правыми частями также не очень сильно отличались друг от друга. К сожалению, так бывает не всегда. Причиной этому - характеристика матрицы, называемая ее обусловленностью . Об этом и пойдет дальше речь.

Сначала требуется ввести меру близости векторов и x, т.е. вектора ошибки . Мерой величины вектора является норма (определяться она может по-разному). Возьмем пока в качестве нормы обычную евклидову норму вектора (квадратный корень из суммы квадратов его компонент), т.е.

Сответственно

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

т.е. спектральная матричная норма равна квадратному корню из максимального собственного числа матрицы A T A. В MATLAB имеется функция norm для вычисления норм матриц и векторов, которая, в частности, умеет вычислять приведенные выше нормы. Почему мы выбрали именно такие нормы векторов и матриц, подробно объясняется в разделе , где приведено немного выкладок и определений. Это связано с оценкой, которой мы будем пользоваться для ошибки в решении системы (вывод этого неравенства также приведен в упомянутом разделе):

Здесь через обозначено число обусловленности матрицы, которое определяется так:

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

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

Для создания матрицы Гильберта в MATLAB предусмотрена функция hilb, входной аргумент которой задает размер матрицы. Возьмем небольшую матрицу 6 на 6:

N = 6; H = hilb(n) H = 1.0000 0.5000 0.3333 0.2500 0.2000 0.1667 0.5000 0.3333 0.2500 0.2000 0.1667 0.1429 0.3333 0.2500 0.2000 0.1667 0.1429 0.1250 0.2500 0.2000 0.1667 0.1429 0.1250 0.1111 0.2000 0.1667 0.1429 0.1250 0.1111 0.1000 0.1667 0.1429 0.1250 0.1111 0.1000 0.0909

X = ones(n, 1); b = H*x b = 2.4500 1.5929 1.2179 0.9956 0.8456 0.7365

Видим, что ни в матрице, ни в векторе правой части нет ничего «подозрительного», все числа не сильно отличаются друг от друга. Теперь сформируем возмущенную правую часть b + δb, добавив в вектор b небольшие числа порядка 10 -5 и решим систему с возмущенной правой частью для получения вектора .

Delta_b = 1e-5*(1:n)"; x_tilda = H\(b + delta_b) x_tilda = 0.9978 1.0735 0.4288 2.6632 -1.0160 1.8593

Видно, что полученное решение далеко от точного, где должны быть все единицы. Вычислим относительную ошибку в решении и в правой части (функция norm по умолчанию вычисляет евклидову норму вектора):

Delta_x = x - x_tilda; LEFT = norm(delta_x)/norm(x) LEFT = 1.1475 RIGHT = norm(delta_b)/norm(b) RIGHT = 2.7231e-005

Итак, ошибка в решении порядка единицы, хотя изменения в правой части были порядка 10 -5 . Это полностью укладывается в приведенное выше неравенство для ошибки. Действительно, вычислим число обусловленности cond(H) при помощи функции MATLAB, которая называется cond и по умолчанию вычисляет для спектральной нормы матрицы

C = cond(H) c = 1.4951e+007

LEFT ≈ 1.1475 меньше

C* RIGHT ? 1.4951e+07 * 2.7231e-05 ≈ 407

и неравенство выполнено (даже с некоторым запасом).

При увеличении размера матрицы Гильберта ошибка в решении будет только расти (это несложно проверить, задавая n = 7, 8, …). Причем, при n = 12 выведется сообщение о том, что матрица плохообусловлена и решение может оказаться неверным

Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 2.409320e-017.

В качестве меры обусловленности здесь выбрана величина RCOND равная единице, деленной на оценку числа обусловленности (число обусловленности оценивается по достаточно быстрому алгоритму, который работает намного быстрее, чем cond, о чем более подробно написано в справке по функции rcond). Если значение RCOND мало, то матрица считается плохообусловленной.

Увеличение ошибки в решении при увеличении размера матрицы Гильберта происходит из-за того, что число обусловленности матрицы Гильберта растет очень быстро с ее размером. В этом несложно убедиться при помощи простого цикла и функции semilogy (масштаб по оси ординат - логарифмический):

N = 1:20; C = zeros(1, 20); for n = N H = hilb(n); C(n) = cond(H); end semilogy(N, C) grid on, title("cond(H)"), xlabel("n")

Кстати, поскольку функция cond находит число обусловленности при помощи численного метода (а именно сингулярного разложения для нахождения сингулярных чисел), то число обусловленности после n = 12 вычисляется уже неверно, на самом деле оно должно расти и дальше, в чем можно убедиться при помощи символьных вычислений в MATLAB и операций с заданным числом значащих цифр

N = 1:20; C = zeros(1, N); digits(60) for n = N H = vpa(sym(hilb(n))); % вычисление матрицы Гильберта с точностью до 60-го знака sigma = svd(H); % нахождение сингулярных чисел матрицы Гильберта % вычисление числа обусловленности матрицы Гильберта C(n) = max(double(sigma))/min(double(sigma)); end semilogy(N, C) grid on, title("cond(H)"), xlabel("n")

Рассмотрим теперь три важных момента, касающихся данного примера.

Первый из них - решение системы Hx = b (с вектором правой части, соответствующим решению из единиц) при помощи оператора обратной косой черты не даст точные единицы, погрешность будет уже в десятом знаке (хотя в MATLAB все вычисления по умолчанию выполняются с двойной точностью)

Format long H\b ans = 0.99999999999916 1.00000000002391 0.99999999983793 1.00000000042209 0.99999999953366 1.00000000018389

Так происходит из-за того, что при вычислении вектора b при умножении матрицы H на вектор из всех единиц мы уже заложили в него некоторую погрешность. Кроме того, ошибки округления в процессе решения системы также сыграли свою роль и плохая обусловленность матрицы (даже небольшого размера) привела к таким ошибкам в решении. Действительно, для матриц небольшого размера с небольшим числом обусловленности такого эффекта наблюдаться не будет. Возьмем матрицу 6 на 6 из случайных чисел, вычислим ее число обусловленности

A = rand(n); cond(A) ans = 57.35245279907571

Затем сформируем правую часть, соответствующую точному решению из всех единиц

X = ones(n, 1); b = A*x;

и решим систему, получив в результате хорошую точность

>> A\b ans = 1.00000000000000 1.00000000000000 1.00000000000000 1.00000000000000 1.00000000000000 1.00000000000000

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

>> det(hilb(6)) ans = 5.3673e-018

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

A = ;

Ее определитель очень мал, порядка 10 -121

Det(A) ans = 9.9970e-121

(слово «очень», конечно условно, но он меньше определителя матрицы Гильберта размера 6 на 6, при решении системы с которой у нас возникли проблемы). Однако, малость определителя никак не повлияет на ошибку в решении при возмущении правой части системы, что несложно показать, сформировав систему с известным решением, внеся возмущения в правую часть и решив систему:

X = ; b = A*x; delta_b = 1e-5*; x_tilda = A\(b+delta_b); delta_x = x - x_tilda; LEFT = norm(delta_x)/norm(x) RIGHT = norm(delta_b)/norm(b) LEFT = 2.1272e-005 RIGHT = 2.1179e-005

Итак, относительная ошибка в решении практически равна относительной ошибке в правой части, поскольку число обусловленности приведенной выше матрицы A немного больше 1, а именно:

C = cond(A) c = 1.0303

И, наконец, рассмотрим третий вопрос, касающийся достижения знака равенства в неравенстве для ошибки в решении

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

в котором относительное возмущение правой части равно 10 -5

RIGHT = norm(delta_b)/norm(b) RIGHT = 1.0000e-005

Относительная ошибка в решении системы оказывается 10 5

X = A\b; x_tilda = A\(b+delta_b); delta_x = x - x_tilda; LEFT = norm(delta_x)/norm(x) LEFT = 1.0000e+005

И происходит так потому, что
1) в этом примере число обусловленности матрицы A равно 10 10 ;

C = cond(A) c = 1.0000e+010

2) неравенство для ошибки в решении переходит в равенство.

Если установить формат long, то мы увидим, что LEFT, RIGHT и число обусловленности не в точности 10 5 , 10 -5 и 10 10 , соответственно, но это вызвано ошибками округлений. При решении в точной арифметике, равенство получилось бы точным, что можно показать аналитически (см. раздел ).

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

|| Ax || V = || b || V ⇒ || A || M || x || V ≥ || b || V

где || || M - некоторая матричная норма. Для того, чтобы это можно было сделать, матричная норма || || M и векторная норма || || V должны удовлетворять следующему неравенству

|| A || M || x || V ≥ || Ax || V

для любых матриц A и векторов x. Если это неравенство выполняется, то матричная норма || || M называется согласованной с векторной нормой || || V . Известно, например, что спектральная норма

(квадратный корень из максимального собственного числа матрицы A T) согласована с евклидовой векторной нормой

Именно поэтому в разделе мы проводили все эксперименты с привлечением данных норм.

Получив неравенство || A || M || x || V ≥ || b || V далее заметим, что из Ax = b и следует, что . Поскольку матрица A невырожденная, то отсюда следует, что δx = A -1 δb и || δx || V = || A -1 δb || V . Опять используем согласованность норм и получаем неравенство || δx || V ≤ || A -1 || M || δb || V . Далее в двух полученных неравенствах

|| A || M || x || V ≥ || b || V и || δx || V ≤ || A -1 || M || δb || V

меньшую величину одного из неравенств делим на большую величину другого неравенства, а большую, соответственно, на меньшую

и простым преобразованием окончательно получаем требуемое неравенство

где cond(A) = || A || M * || A -1 || M .

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

согласована с первой векторной нормой

а максимальная строчная матричная норма, вычисляемая по формуле

согласована с бесконечной векторной нормой

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

  • q = norm(A, 1) - максимальная столбцевая норма матрицы A
  • q = norm(A, inf) - максимальная строчная норма матрицы A
  • q = norm(a, 1) - первая векторная норма a
  • q = norm(a, inf) - бесконечная векторная норма a

Число обусловленности матрицы cond(A) = || A || M * || A -1 || M относительно различных матричных норм может быть вычислено при помощи функции cond. Если cond вызывается с одним входным аргументом (матрицей), то вычисляется число обусловленности относительно спектральной матричной нормы. При вызове функции cond с дополнительным аргументом возвращаются числа обусловленности относительно указанной матричной нормы:

  • с = cond(A, 1) - число обусловленности относительно максимальной столбцевой нормы матрицы
  • с = cond(A, inf) - число обусловленности относительно максимальной строчной нормы матрицы

В качестве примера, демонстрирующего важность согласованности норы матрицы с нормой вектора в неравенстве для ошибки, приведем пример с матрицей A, вектором правой части b и вектором ошибок в правой части delta_b.

A = ; b = [ -5.7373057243726833e-001 -1.5400413072907607e-001 -5.3347697688693385e-001 -6.0209490373259589e-001]; delta_b = [-0.71685839091451e-5 0.54786619630709e-5 0.37746931527138e-5 0.20850322383081e-5];

Вычислим правую и левую часть оценки с использованием первой векторной нормы, а число обусловленности матрицы по отношению к спектральной матричной норме

RIGHT = norm(delta_b,1)/norm(b,1); c = cond(A); x = A\b; x_tilda = A\(b+delta_b); delta_x = x - x_tilda; LEFT = norm(delta_x,1)/norm(x,1); format short e disp()

Получаем для левой и правой частей неравенства следующие значения
9.9484e+004 9.9323e+004

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

Разберем теперь, почему число обусловленности матрицы не может быть меньше единицы. Число обусловленности матрицы A определяется как cond(A) = || A || M * || A -1 || M , где || A || M - какая-то матричная норма A. Матричная норма (т.е. правило, по которому каждой матрице сопоставляется число) не может быть произвольной, она должна удовлетворять следующим четырем аксиомам:

A1. || A || ≥ 0 для любой матрицы A и || A || = 0 тогда и только тогда, когда A - нулевая матрица.

А2. || αA || = | α | * || A || для любой матрицы A и числа α.

А3. || A + B || ≤ || A || + || B || для любых матриц A и B

А4. || AB || ≤ || A || * || B || для любых матриц A и B.

Из последней аксиомы видно, что норма определена только для квадратных матриц (хотя, в приведенных выше формулах для вычисления различных норм, в принципе, нет такого ограничения). Кроме того, из последней аксиомы следует, что любая норма единичной матрицы I не меньше единицы, действительно

|| I || = || I*I || ≤ || I || 2 ⇒ || I || ≥ 1.

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

1 ≤ || I || = || AA -1 || ≤ || A || * || A -1 || = cond(A).

Завершим данный раздел исследованием причины появления точного равенства в оценке ошибки в решении при возмущении правой части:

и построением соответствующих примеров в MATLAB. (Приводимое ниже рассуждение содержится, например, в книге Дж. Форсайт, К Молер. Численное решение систем линейных алгебраических уравнений. М: «Мир», 1969.)

Существенную роль в рассуждениях играет теорема о сингулярном разложении матрицы, согласно которой для любой вещественной матрицы размера на существуют две ортогональные матрицы U и V размера n на n (U T U=UU T и V T V = VV T) такие, что произведение D = U T AV является диагональной матрицей, причем можно выбрать U и V так, что

где μ 1 ≥ μ 2 ≥…≥μ r ≥ μ r+1 =…=μ n =0,

а r - ранг матрицы A. Числа μ k называют спектральными числами матрицы A. Для невырожденной матрицы A верно:

μ 1 ≥ μ 2 ≥ … ≥μ n ≥ 0.

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

|| Ux || = || x ||.

Поскольку умножение на ортогональную матрицу не изменяет спектральной нормы, то

следовательно для числа обусловленности матрицы верно равенство

У нас есть две системы: Ax = b (с точным решением x) и (с точным решением ). Очевидно, что ошибка δx является решением системы, правая часть которой является возмущением δb, т.е. системы Aδx = δb. Пусть D = U T AV есть сингулярное разложение матрицы A, тогда UDV T = A из-за того, что U и V - ортогональные матрицы. Далее

Ax = b ⇒ UDV T x = b.

Введем обозначения

x" = V T x, b" = U T b,

тогда следующие системы эквивалентны

Ax = b ⇔ Dx" = b"

Совершенно аналогично рассмотрим систему относительно ошибки

Aδx = δb ⇒ UDV T δx = δb

Вводим обозначения

δx" = V T δx, δb" = U T δb,

при которых следующие системы оказываются эквивалентны

Aδx = δb ⇔ Dδx" = δb"

Т.е. мы получили простые эквивалентные системы с диагональной матрицей, на диагонали которой стоят спектральные числа матрицы A.

Выберем теперь специальным образом правые части этих систем, именно, пусть

где β > 0, тогда соответствующее ей решение системы Dx" = b" будет

где μ 1 - максимальное сингулярное число матрицы A. Схожим образом поступим и с системой Dδx" = δ b" , именно, пусть

где γ > 0, тогда соответствующее ей решение системы Dδx" = δb" будет

Из сохранения нормы вектора при его перемножении с ортогональной матрицей следует, что

β/μ 1 = || x" || = || V T x || = || x || и γ/μ n = || δx" || = || V T δx || = ||δx ||.

Совершенно аналогично получаем равенства:

β = || b" || = || U T b || = || b || и γ = || δb" || = || U T δb || = || δb ||.

а поскольку

то окончательно получаем

Итак, для каждой матрицы A можно построить вектор правой части системы и его возмущение такие, что ошибка в решении будет произведением числа обусловленности матрицы на ошибку в правой части. В MATLAB соответствующее построение может быть выполнено и без использования сингулярного разложения (хотя оно может быть найдено при помощи функции svd).

Сначала задаем n и получаем две ортогональные матрицы U и V, выполняя QR-разложение матриц n на n из случайных чисел:

N = 4; = qr(rand(n)); = qr(rand(n));

D = diag();

После этого конструируем матрицу A, используя диагональную матрицу D и ортогональные матрицы U и V

A = U*D*V";

Число обусловленности матрицы A будет совпадать с числом обусловленности матрицы D, которое в нашем примере равно 10 10

Beta = 1; gamma = 1e-5;

и строим вектора b" и δb"

B1 = "; db1 = ";

по которым находим вектора b и δb

X = A\b; x_tilda = A\(b+db);

вычисляем левую и правую части неравенства

dx = x - x_tilda; RIGHT = norm(db)/norm(b); LEFT = norm(dx)/norm(x);

и выводим их

Format short e disp()

Получаем равенство

Министерство образования Республики Беларусь

Учреждение образования

· contour 3(X , Y , Z ) – строит контурные линии для поверхности, полученные при расслоении трехмерной фигуры рядом секущих плоскостей, расположенных параллельно опорной плоскости фигуры.

В приведенном ниже примере показаны возможности применения описанных команд для построения графика поверхности Z = sin(X) cos(X).

>> =meshgrid(-3:0.1:3,-3:0.1:3); Z=sin(X).*cos(X);

>> subplot(3,2,1), plot3(X, Y,Z) % Рисунок 4.3 (а)

>> subplot(3,2,2), mesh(X, Y,Z) % Рисунок 4.3 (б)

>> subplot(3,2,3), surf(X, Y,Z) % Рисунок 4.3 (в)

>> subplot(3,2,4), surfc(X, Y,Z) % Рисунок 4.3 (г)

>> subplot(3,2,5),meshz(X, Y,Z) % Рисунок 4.3 (д)

>> subplot(3,2,6),contour3(X, Y,Z) % Рисунок 4.3 (е)


4.3. Форматирование графиков

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

Таблица 4.2 - Форматирование графиков

Действие

Средства графического интерфейса

Переход в режим редактирования.

Щелкнуть по кнопке Edit Plot в панели инструментов графического окна.

Форматирование линий и маркеров опорных точек графиков.

В режиме редактирования дважды щелкнуть по линии графика левой клавишей мыши. В появившемся окне Property Editor– Line установить все необходимые параметры линии (толщину, стиль, цвет и т. д.).

plot(X, Y,S), plot3(X, Y,Z, S)

(Описание команд приводится в п.4.1 и п.4.2)

Форматирование осей графиков.

В режиме редактирования дважды щелкнуть по оси графика. В появившемся окне Property Editor– Axes установить все необходимые параметры осей.

Заголовок графика и метки осей можно также с помощью команд Insert Title , Insert Xlabel, Insert Ylabel главного меню графического окна.

axes – управляет свойствами осей.

grid - включает и отключает координатную сетку.

xlabel(S), ylabel(S), zlabel(S) – устанавливает надписи возле осей. Здесь S – строковая константа или переменная.

title(S) – выводит заголовок графика

Нанесение надписей прямо на график.

Щелкнуть по кнопке Insert Text ,зафиксировать место надписи щелчком мыши и ввести нужный текст.

text(X, Y, S) – добавляет в двумерный график текст, заданный строкой S, так что начало текста расположено в точке с координатами (X, Y).

text(X, Y, Z, S) - добавляет текст в трехмерный график.

Нанесение линий и линий со стрелками прямо на график

Щелкнуть по одной из кнопок Insert Arraw или Insert Line . Установить указатель мыши в нужное место графика и при нажатой левой клавише мыши нарисовать линию.

Построение легенды

Insert, а затем команду legend.

legend(S1, S2, S3,…) – добавляет к текущему графику легенду с пояснениями в виде строк, указанных в списке параметров.

Вывод шкалы цветов

В главном меню графического окна выбрать команду Insert, а затем команду Colorbar.

colorbar(‘ vert’), colorbar(‘ horiz’) – выводит вертикальную или горизонтальную шкалу цветов.

Вращение графика

Щелкнуть по кнопке Rotate 3 D и вращать график с помощью мыши (может применяться также для двумерных графиков).

rotate3 d – задает поворот трехмерной фигуры.

с граничными условиями y (t 0 , t end , p ) = y , где t end , t 0 начальные и конечные точки интервалов. Параметр t (независимая переменная) необязательно означает время, хотя чаще всего решение ДУ ищется во временной области. Система ДУ в форме Коши записывается аналогично (1.1), но под y в этом случае подразумевается вектор-столбец зависимых переменных. Вектор p задает начальные условия.

Для решения ДУ второго и высшего порядка их нужно свести к системе ДУ первого порядка.

Возможны ДУ, не разрешенные относительно производной:

F (t , y , dy /dt ) = 0. (1.2)

Уравнения (1.2) аналитически к форме (1.1) обычно привести не удается. Однако численное решение особых трудностей не вызывает достаточно для определения f (y , t ) решить (1.2) численно относительно производной при заданных y и t .

Решатели ОДУ

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

В этом разделе обобщенное название solver (решатель) означает один из возможных численных методов решения ОДУ: ode45, ode23, ode113, ode15s, ode23s, ode23t, ode23tb, bvp4c или pdepe.

Решатели реализуют следующие методы решения систем ДУ:

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

Ode23 одношаговые явные методы Рунге-Кутта 2-го и 4-го порядков в модификации Богацки и Шампина. При умеренной жесткости системы ОДУ и низких требованиях к точности этот метод может дать выигрыш в скорости решения.

Ode113 многошаговый метод Адамса-Башворта-Мултона переменного порядка класса предиктор-корректор. Это адаптивный метод, который может обеспечить высокую точность решения.

Ode15s многошаговый метод переменного порядка (от 1 до 5, по умолчанию 5), использующий формулы численного «дифференцирования назад». Это адаптивный метод, его стоит применять, если решатель ode45 не обеспечивает решения и система ДУ жесткая.

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

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

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

Bvp4c служит для проблемы граничных значений систем ДУ вида y ′ = f (t , y ), F (y (a ), y (b ), p ) = 0 (полная форма системы уравнений Коши). Решаемые им задачи называют двухточечными краевыми задачами, поскольку решение ищется при задании граничных условий как в начале, так и в конце интервала решения.

Все решатели могут решать системы уравнений явного вида y ′ = F (t , y ), причем для решения жестких систем уравнений рекомендуется использовать только специальные решатели ode15s, ode23s, ode23t, ode23tb.

Использование решателей систем ОДУ

tspan вектор, определяющий интервал интегрирования [t 0 t final ]. Для получения решений в конкретные моменты времени t 0 , t 1 , …, t final (расположенные в порядке уменьшения или увеличения) нужно использовать tspan = [ t 0 t 1 … t final ];

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

Options аргумент, создаваемый функцией odeset (еще одна функция odeget или bvpget (только для bvp4c) позволяет вывести параметры, установленные по умолчанию или с помощью функции odeset/bvpset);

p 1, p 2,… произвольные параметры, передаваемые в функцию F ;

T , Y матрица решений Y , где каждая строка соответствует времени, возвращенном в векторе-столбце T .

Перейдем к описанию синтаксиса функций для решения систем ДУ (под именем solver подразумевается любая из представленных выше функций).

[T ,Y ]=solver(@F ,tspan ,y 0) интегрирует систему ДУ вида y ′ = F (t , y ) на интервале tspan с начальными условиями y 0 . @F дескриптор ОДУ-функции (можно также задавать функцию в виде "F "). Каждая строка в массиве решений Y соответствует значению времени, возвращаемому в векторе-столбце T .

[T ,Y ]=solver(@F ,tspan ,y 0 ,options) дает решение, подобное описанному выше, но с параметрами, определяемыми значениями аргумента options, созданного функцией odeset. Обычно используемые параметры включают допустимое значение относительной погрешности RelTol (по умолчанию 1e3) и вектор допустимых значений абсолютной погрешности AbsTol (все компоненты по умолчанию равны 1e6).

[T ,Y ]=solver(@F ,tspan ,y 0 ,options,p 1 ,p 2 …) дает решение, подобное описанному выше, передавая дополнительные параметры p 1 , p 2 , … в m -файл F всякий раз, когда он вызывается. Используйте options=, если никакие параметры не задаются.

Решение ОДУ первого порядка

ПОРЯДОК ВЫПОЛНЕНИЯ РАБОТЫ

· титульный лист;

· исходные данные варианта;

· решение задачи;

· результаты решения задачи.

Пример

Найти решение дифференциального уравнения на отрезке , для которого у (1,7) = 5,3.

Создаем в Command Window функция пользователя

g=@(x,y);

В синтаксисе функции @(x,y) x независимая переменная, y зависимая переменная, x -cos(y /pi ) правая часть ДУ.

Процесс решения осуществляется обращением в Command Window к решателю (солверу) следующим оператором:

Ode23(g,,);

Построение графика с сеткой осуществляется следующими операторами:

Результат представлен на рис. 1.1

Рис. 1.2.1. Визуализация численного решения

ЗАДАНИЕ

1. Найдите решения ДУ первого порядка , удовлетворяющего начальным условиям у(х 0 ) = у 0 на промежутке [a , b ].

2. Построить графики функции.

Варианты заданий .

№ варианта у(х 0 )=у 0 [a , b ]
y 0 (1,8)=2,6
y 0 (0,6)=0,8
y 0 (2,1)=2,5
y 0 (0,5)=0,6
y 0 (1,4)=2,2
y 0 (1,7)=5,3
y 0 (1,4)=2,5
y 0 (1,6)=4,6
y 0 (1,8)=2,6
y 0 (1,7)=5,3
y 0 (0,4)=0,8
y 0 (1,2)=1,4

Лабораторная работа № 2

Решение систем ОДУ

ЦЕЛЬ РАБОТЫ

Сформировать у студентов представления о применении систем ДУ в различных областях; привить умения решать задачу Коши для систем ДУ.

ПОРЯДОК ВЫПОЛНЕНИЯ РАБОТЫ

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

2. Оформите отчет по лабораторной работе, который должен содержать:

· титульный лист;

· исходные данные варианта;

· решение задачи;

· результаты решения задачи.

Пример

Решить систему

с использованием решателя ode23().

Решение:

1. Создать в редакторе m-файл функции вычисления правых частей ДУ.

Пусть имя в редакторе файла sisdu.m, тогда функция может иметь следующий вид:

function z=sisdu(t,y)

z1=-3*y(2)+cos(t)-exp(t);

z2=4*y(2)-cos(t)+2*exp(t);

>> t0=0;tf=5;y0=[-3/17,4/17];

>> =ode23("sisdu",,y0);

>> plot(t,y)

Рис. 1.3.1. Визуализация численного решения, полученного с помощью функции ode23.

1. Что значит решить задачу Коши для системы ДУ?

2. Какие существуют методы решения систем ДУ?

ЗАДАНИЕ

1. Найдите решение системы ДУ

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

2. Построить графики функций.

Для примера приводится функция решения 8-го варианта:

function z=ssisdu(t,y)

% вариант 8

z1=-a*y(1)+a*y(2);

z2=a*y(1)-(a-m)*y(2)+2*m*y(3);

z3=a*y(2)-(a-m)*y(3)+3*m*y(4);

z4=a*y(3)-3*m*y(4);

>> =ode23("ssisdu",,);

>> plot(t,100*y)

Рис. 1.3.2. Визуализация численного решения, полученного с помощью функции ode23.

Варианты заданий .

№ варианта Задания
a m
0,1 1,2
0,2 1,5
0,3 1,7
0,4 1,9
0,5
0,6 1,9
0,7 2,3
0,8 2,7
0,9
0,1 1,5
0,2 1,1
0,3

Лабораторная работа № 3

1.4Решение ОДУ n -го порядка

ЦЕЛЬ РАБОТЫ

Сформировать у студентов представления о применении ДУ высших порядков в различных областях; привить умения решать задачу Коши для ДУ высших порядков с помощью прикладных программ; развить навыки проверки полученных результатов.

ПОРЯДОК ВЫПОЛНЕНИЯ РАБОТЫ

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

2. Оформите отчет по лабораторной работе, который должен содержать:

· титульный лист;

· исходные данные варианта;

· решение задачи;

· результаты решения задачи.

Пример 1.

Решить ДУ второго порядка при данных начальных условиях .

Решение:

Сначала приведем ДУ к системе:

1. Создать m-файл функции вычисления правых частей ДУ.

Пусть имя файла sisdu_3.m, тогда функция может иметь следующий вид:

function z=sisdu_3(x,y)

z2=6*x*exp(x)+2*y(2)+y(1);

2. Выполнить следующие действия:

>> x0=0;xf=10;y0=;

>> =ode23("sisdu_3",,y0);

>> plot(x,y(:,1))

Рис. 1.4.1. Визуализация численного решения, полученного с помощью функции ode23.

ПРИМЕРНЫЕ ВОПРОСЫ НА ЗАЩИТЕ РАБОТЫ

1. Что значит решить задачу Коши для ДУ высших порядков?

2. Как привести ДУ m -го порядка к системе ДУ?

ЗАДАНИЕ

1. Найдите решение ДУ, удовлетворяющее начальным условиям на промежутке .

2. Построить графики функций.

Варианты заданий .

№ варианта Задания
Уравнения Начальные условия







Лабораторная работа № 4 – 5

Динамические системы (ДС)

ЦЕЛЬ РАБОТЫ

Знакомство студентов с основными понятиями ДС, их классификация, фазовое пространство ДС, кинематическая интерпретация системы ДУ, эволюция ДС. Уравнение движения маятника. Динамика осциллятора Ван дер Поля.

2. Динамическая система (ДС) математический объект, соответствующий реальным системам (физическим, химическим, биологическим и др.), эволюция которых однозначно определяется начальным состоянием. ДС определяется системой уравнений (дифференциальных, разностных, интегральных и т.д.), допускающих существование на бесконечном интервале времени единственность решения для каждого начального условия.

Состояние ДС описывают набором переменных, выбираемых из соображений естественности их интерпретации, простоты описания, симметрии и т.п. Множество состояний ДС образует фазовое пространство, каждому состоянию отвечает точка в нём, а эволюция изображается (фазовыми) траекториями. Чтобы определить близость состояний, в фазовом пространстве ДС вводят понятие расстояния. Совокупность состояний в фиксированный момент времени характеризуется фазовым объёмом.

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

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

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







2024 © gtavrl.ru.