Решение систем дифференциальных уравнений

Система MLимеет встроенные средства для решения обыкновенных дифференциальных уравнений (ДУ) любого порядка или систем ДУ. Дифференциальным уравнением n-го порядка называется соотношение вида F(t, y, y’,…, y⁽ⁿ⁾)=0. Наивысший порядок производной искомой функции в данном уравнении называется порядком дифференциального уравнения. Решением ДУ является такая функция, которая при подстановке в него обращает его в равенство. В случае системы дифференциальных уравнений n-го порядка решением являются n функций. При использовании численных методов решения искомые функции выражаются не аналитически, а в табличном виде (в виде набора значений в определенных точках). Самое простое ДУ – это уравнение 1-го порядка, разрешенное относительно производной. Оно имеет вид y’=f(t, y). ДУ и системы 2-го порядка и выше можно свести к системам 1-го порядка. Так, для того чтобы свести к системе 1-го порядка уравнение y’’=f(t, y, y’), нужно использовать замену: y₁=y и y₂=y’. В результате замены будем иметь систему

В ML имеется целый ряд встроенных функций, предназначенных для решения систем дифференциальных уравнений. Это функции ode23, ode45, ode113. В них реализованы различные методы решения обыкновенных ДУ. Наиболее часто используют функции ode23 и ode45, в основе работы которых лежит метод Рунге–Кутта. Для систем уравнений 2-го и 3-го порядков используют функцию ode23, а для систем 4-го и 5-го порядков – функцию ode45. Обращение к функции, находящей решение системы ДУ, имеет вид

[t, y]= <имя встроенной функции>(fun, interval, y0 ),

где fun – имя m-файла, в котором вычисляются правые части системы ДУ; interval – вектор, определяющий отрезок, на котором ищется решение; y0 – скаляр или вектор-столбец, в котором задаются начальные условия; t – вектор-столбец, содержащий значения независимой переменной на заданном интервале (обычно это время t); y – матрица решений. Количество столбцов в ней равно порядку системы ДУ. Каждый столбец содержит значения одной искомой функции на заданном отрезке. Каждая строка соответствует определенному значению t.

П р и м е р. Имеем ДУ 2-го порядка y’’+4y’+5y=0, описывающее движение точки в зависимости от времени. Требуется его решить на отрезке [0,5] с начальными условиями y(0)=3 и y’(0)= -5.

Выполним замену: y₁=y и y₂=y’, получим систему ДУ:

при y₁(0)=3 и y₂(0)=-5 (тогда y₁ будет обозначать координаты точки, а y₂ – ее скорость).

Составим файл-функцию для вычисления правых частей системы ДУ. Эта функция должна иметь два входных параметра: переменную t, вектор y с количеством элементов, равным числу неизвестных системы, и один выходной параметр (вектор правых частей системы):

function F=mydif(t,y)

F=[y(2); -5*y(1)-4*y(2)];

Сохраним файл с именем mydif в каталоге, разрешенном для записи. Обратимся к функции ode45:

>>[t,y]=ode45('mydif', [0 5],[3;-5])

В результате будут получены вектор-столбец t – вектор значений аргумента (в нашем примере это время) и матрица y из двух столбцов. Первый столбец содержит значения функции y₁, т.е. самой искомой функции при каждом значении t, а второй столбец – значения функции y₂, т.е. первой производной искомой функции.
У нас в первом столбце будут получаться значения координат точки, а во втором – значение скорости ее движения в каждый момент времени t.

Отобразим график решения ДУ и график производной (рис. 34):

>> plot(t,y(:, 1), 'r', t, y (:, 2),'k--')

>> xlabel('t');

>> ylabel('y')

Рис. 34

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

Знакомство с пакетом MAtlab,
решение простейших задач в системе MAtlab

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

Постановка задачи

1. Вычисление значений функции на заданном отрезке.

2. Выполнение действий с векторами.

3. Формирование и выполнение действия с матрицами.

4. Построение графиков функций одной переменной на заданном отрезке.

5. Построение графика функции двух переменных.

6. Создание скрипт-файла.

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

Индивидуальные задания

Задание №1.Вычислить N значений функции на заданном отрезке. На экран вывести значения аргумента и значения функции.

Варианты

Функция Отрезок Количество разбиений

1. [0,2 ] N=10

2. [-0.2,4] N=9

3. [0,0.3] N=7

4. [0,1] N=10

5. [0,3] N=8

6. [ ,3 ] N=8

7. [-1,1] N=7

8. [-1,1] N=10

9. [-2,2] N=7

10. [-2,2] N=9

Задание №2.Для заданных векторов a и b длиной n (значения элементов векторов и их длину студент задает сам) выполнить преобразования и вычисления в соответствии с вариантом.

Варианты

1. В векторе a элементы с номерами от n1 до n2 удвоить, а в векторе b заменить их средним арифметическим.

2. Образовать новый вектор c=[a₁,a₂,…,a_n,b₁,b₂,…,b_n], определить его максимальный и минимальный элементы и поменять их местами.

3. Образовать вектор c =[a₁,a₂,a₃,b₄,b₅,…,b_n] и упорядочить его по возрастанию и убыванию.

4. Образовать вектор c =[a₃,a₄,…,a_n,b₁,b₂,b₃] и переставить элементы вектора c в обратном порядке. Результат записать в новый вектор.

5. Получить вектор x, содержащий удвоенные значения элементов вектора a, и вектор y, содержащий утроенные значения элементов вектора b. Определить среднее арифметическое каждого вектора.

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

7. Получить два новых вектора, состоящих из элементов исходных векторов, начиная с номера n1 до номера n2. Найти сумму минимальных элементов новых векторов.

8. Заменить нулем минимальный элемент вектора a и максимальный элемент вектора b.

9. Вычислить произведение элементов векторов с номерами от n1 до n2. Найти минимальные значения векторов и заменить последние элементы векторов их минимумами.

10. Образовать вектор c =[a₂,a₃,a₄,b₃,b₄,…,b_n]. Элементы с номерами от n1 до n2 заменить средним арифметическим этих элементов.

Задание №3.При помощи встроенных функций для заполнения стандартных матриц, индексации двоеточием и, возможно, объединения, поворота или транспонирования получить следующие матрицы. Применить функции обработки данных и поэлементные операции для нахождения заданных величин.

Варианты

1. A=

2. A= .

3. A= .

4. A= .

5. A= .

6. A= .

7. A= .

8. A= .

9. A= .

10. A= .

Задание №4.Построить графики двух функций на заданном отрезке. Вывести графики:

• в разных окнах,

• в одном окне в одних осях,

• в одном окне в разных осях.

Использовать различные цвета, стили, подписи, легенду. Нанести сетку.

Варианты

Функция f Функция g Аргумент x

1.

2.

3.

4.

5.

6.

7.

8.

9.

10.

Построить график кусочно-заданной функции, отобразить ветви разными цветами и маркерами.

Варианты

1

2

3

4

5

6

7

8

9

10

Задание №5.Построить график функции двух переменных.

Варианты

Функция z Аргумент x Аргумент y

1.

2.

3.

4.

5.

6.

7.

8.

9.

10

Задание №6.Написать скрипт-файл для решения следующих задач.

Варианты

1. По заданному вектору определить номер его элемента с наибольшим отклонением от среднего арифметического всех элементов векторов.

2. Найти среднее арифметическое элементов заданного вектора и заменить первый элемент этим значением.

3. Вычислить максимальное значение среди элементов главной диагонали заданной матрицы.

4. Переставить первый столбец квадратной матрицы со строкой, где находится наименьший элемент.

5. Сложить все элементы заданной матрицы, кроме элементов главной диагонали.

6. Заменить максимальный элемент вектора средним значением всех его элементов.

7. Заменить элемент матрицы с индексами 1,1 произведением всех элементов матрицы.

8. Заменить последний элемент вектора максимальным элементом.

9. Найти значение и номер максимального элемента. Графически отобразить элементы заданного вектора синими маркерами, а максимальный элемент – красным.

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

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

Решение типовых вычислительных
задач в системе MAtlab

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

Постановка задачи

1. Вычисление корней полинома.

2. Решение системы линейных уравнений.

3. Вычисление локального минимума и максимума функции.

4. Вычисление определенного интеграла.

5. Решение трансцендентного уравнения.

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

· для каждого пункта задания создать свой скрипт-файл;

· написать файл-сценарий (script-файл), в котором пользовательский интерфейс оформлен в виде меню. Выбранный пункт меню определяет выполнение файла, соответствующего пункту индивидуального задания. При написании программы для реализации меню использовать встроенную функцию menu(), оператор цикла с предусловием while и функцию eval();

· уравнения и функции для выполнения пунктов № 3–5 описать в виде файлов-функций;

· для пункта №1 при вычислении корней полинома построить график полинома и отобразить на нем найденные действительные корни. Ввод границ построения графика должен осуществляться с клавиатуры;

· для пункта №2 при решении системы линейных уравнений осуществить проверку полученного решения;

· для пункта №3 при поиске максимума и минимума функции построить график заданной функции в заданных границах и отобразить на нем полученные экстремумы маркерами разного цвета;

· для пункта №4 при вычислении интеграла построить график подынтегральной функции, границы графика вводить с клавиатуры и закрасить площадь, ограниченную функцией на заданном отрезке;

· для пункта №5 при решении трансцендентного уравнения построить график функции в границах, заданных пользователем в форме диалога, и отобразить на нем значение корня уравнения цветом, отличным от цвета графика;

· на всех графиках должны быть выведены заголовки, названия осей координат. Все графики должны располагаться в одном графическом окне.

При написании программ обязательно:

· использовать комментарии, содержащие назначение программы и описание ее переменных;

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

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

Варианты индивидуальных заданий

Вариант 1

1. Вычислить корни полинома

2. Найти решение системы линейных уравнений

3. Найти значение локального минимума и максимума функции

4. Вычислить значение определенного интеграла

5. Решить трансцендентное уравнение

Вариант 2

1. Вычислить корни полинома

2. Найти решение системы линейных уравнений

3.

4. Вычислить значение определенного интеграла

5. Решить трансцендентное уравнение

Вариант 3

1. Вычислить корни полинома.

2. Найти решение системы линейных уравнений

3. Найти значение локального минимума и максимума функции

4. Вычислить значение определенного интеграла

5. Решить трансцендентное уравнение

Вариант 4

1. Вычислить корни полинома.

2. Найти решение системы линейных уравнений

3. Найти значение локального минимума и максимума функции

4. Вычислить значение определенного интеграла

5. Решить трансцендентное уравнение

Вариант 5

1. Вычислить корни полинома

2. Найти решение системы линейных уравнений

3. Найти значение локального минимума и максимума функции

4. Вычислить значение определенного интеграла

5. Решить трансцендентное уравнение

Вариант 6

1. Найти корни полинома.

2. Найти решение системы линейных уравнений

3. Найти значение локального минимума и максимума функции

4. Вычислить значение определенного интеграла

5. Решить трансцендентное уравнение

Вариант 7

1. Вычислить корни полинома

2. Найти решение системы линейных уравнений

3. Найти значение локального минимума и максимума функции

4. Вычислить значение определенного интеграла

5. Решить трансцендентное уравнение

Вариант 8

1. Вычислить корни полинома

2. Найти решение системы линейных уравнений

3. Найти значение локального минимума и максимума функции

4. Вычислить значение определенного интеграла

5. Решить трансцендентное уравнение

Вариант 9

1. Вычислить корни полинома

2. Найти решение системы линейных уравнений

3. Найти значение локального минимума и максимума функции

4. Вычислить значение определенного интеграла

5. Решить трансцендентное уравнение

Вариант 10

1. Вычислить корни полинома

2. Найти решение системы линейных уравнений

3. Найти значение локального минимума и максимума функции

4. Вычислить значение определенного интеграла

5. Решить трансцендентное уравнение

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

Программирование в системе MAtlab

Цель работы – освоить операторы языка ML, научиться использовать их для решения ситуационных задач.

Постановка задачи

Написать скрипт-файл для решения задачи из области математики или физики.

Варианты индивидуальных заданий

Вариант 1. В гелиоцентрической системе отсчета Земля движется по окружности радиуса R₁=1,496·10⁸ км(период обращения Т₁=3,156·10⁷ с). Координаты Земли описываются зависимостями

.

Луна, в свою очередь движется вокруг Земли по окружности радиусом R₂=3,844·10⁵ км (период обращения T₂=2,36·10⁶ c). Координаты Луны в геоцентрической системе координат:

.

Построить орбиты Земли и Луны в гелиоцентрической системе координат. Промоделировать, как изменится картина при других значениях R₂ и T₂, например при R₂=3,844·10⁷ км, Т₂=2,36·10⁵с.

Вариант 2. Получить амплитудно-частотную и фазочастот-ную характеристики (АЧХ и ФЧХ) цифрового рекурсивного фильтра N-го порядка:

;

,

где ; ; ; ; и – весовые коэффициенты фильтра,
– период дискретизации. Принять = 0,001 с, N = 2,
– изменять от 0 до π.

Значения весовых коэффициентов вводить с клавиатуры:

a₀=1; a₁ = -2,208; a₂ = 1,208; b₁ = -0,848; b₂ = 0,36.

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

Математическая модель закона изменения ускорения:

Получить выражения для законов изменения скорости v(t) и пути S(t) и построить зависимости g(t), v(t), S(t) для следующих исходных данных:

.

Вариант 4. Построить траекторию спуска космического аппарата в трехмерном пространстве.

Законы изменения составляющих ускорения:

Путем двойного интегрирования получить законы изменения координат x(t), y(t)и h(t) для следующих начальных условий:

Исходные данные:

g_xm = g_ym = g; g_hm = -4g; g = 9,8 м/с²; t_k = 5 мин; шаг изменения времени ∆t = 3 c.

Воспользоваться оператором plot3, местоположение указывать зеленой звездочкой. С помощью операторов line показать проекции точек положения космического аппарата на горизонтальную и вертикальную плоскости.

Вариант 5. Построить траекторию стартующей ракеты в трехмерном пространстве. Законы изменения координат:

при 0 ≤ t ≤10;

при 10 < t ≤ 50,

где

Воспользоваться оператором plot3, местоположение ракеты указывать красной звездочкой. С помощью операторов line показать проекции точек положения летательного аппарата на горизонтальную и вертикальную плоскости. Шаг изменения времени ∆t = 1 c.

Вариант 6. Изобразить интерференционную картину, получившуюся при освещении оранжевым светом с длиной волны
l = 0,6 мкм плоской пластины с прижатыми к ней плосковыпуклыми линзами с радиусом кривизны выпуклой поверхности R = 5м.

Разность фаз интерферирующих волн: , где r – расстояние до центра интерференционной картины. r = 0÷0,4·10^{-2 м (шаг изменения} Dr = 0,5·10^-4 м). Яркость .

Построить график зависимости I(r) в декартовой и полярной системах координат. Для моделирования интерференционной картины ограничить м и построить матрицу значений фазы и интенсивности отраженного света размером 60×60.

; ,

где координаты и изменяются так:

, , N=60.

Построить график зависимости интенсивности отраженного света от координат и карту линий одного уравнения (командой Contour).

Вариант 7. Биоритмы человека представляют собой синусоиды, выходящие из нуля в день рождения человека и имеющие периоды: интеллектуальный – 33 дня, эмоциональный – 28 дней, физический – 23 дня.

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

Вариант 8. Колесо электровоза, движущегося со скоростью , имеет радиус R = 1 м. Необходимо рассчитать и построить траекторию точки, лежащей на расстоянии r = 0,5 м от оси колеса. Считать, что в начальный момент времени точка находилась в самом нижнем положении. Кинематические уравнения движения точки:

, .

Построить график на интервале . Указать на графике положения точки в моменты и . Изобразить вектора скорости движения точки для моментов и .

Вариант 9. Осуществить гармонический синтез пилообразного сигнала по первым 3, 6и 15гармоникам:

.

Для этого суммировать 3, 6 и 15 синусоидальных сигналов
соответственно. Построить графики полученных сигналов при
T= 50, t=0 ÷ T.

Вариант 10.В момент преодоления самолетом звукового барьера число Маха становится равным единице. Функция, определяющая число Маха:

,

где v – скорость полета самолета; H – высота полета; коэффициенты: a = 1222,5; b = 6,875·10^-6; c = 0,3048; d = -5,2656; e = 0,286.

Получить графики зависимости M(v) для H=500;1000; 2000;5000; показать уровень М=1, соответствующий достижению скорости звука.

Из условия M(v,H)=1 получить зависимость скорости преодоления звукового барьера v от высоты H. Для этого, изменяя H в диапазоне от 0 до 2,5·10⁴, решать уравнение M(v,H) – 1 = 0. При решении уравнения передавать H в функцию, описывающую правую часть уравнения, как глобальное данное (командой global H). Построить график зависимости v(H), при котором M(v,H) = 1.

(Для сведения: при ; при ).

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

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

1. Каково назначение системы Matlab?

2. В каких режимах может выполняться работа в системе ML?

3. В виде каких файлов хранится большинство команд и функций системы ML?

4. Опишите структуру окна рабочей среды ML.

5. Для чего в ML используют клавиши управления ­ и ¯?

6. Что является элементарной единицей данных языка ML?

7. Как записываются действительные числа в ML ?

8. Какими командами можно получить информацию о данных, хранящихся в рабочем пространстве?

9. Какие форматы вывода числовых данных в ML вы знаете?

10. Как изменить формат вывода результатов вычисления в ML?

11. Как в системе MLопределяется тип переменных?

12. Назовите правила составления имен переменных.

13. Какие основные системные переменные существуют в ML?

14. Какие операции существуют вML?

15. Назовите операции ML в порядке убывания приоритета.

16. Как представляются вектора и матрицы в ML?

17. Как записываются и чем различаются матричные и поэлементные операции в ML?

18. Назовите способы задания векторов в ML.

19. Как обратиться к элементам векторов и матриц?

20. Для чего используется функция length()?

21. Как обратиться к последнему элементу вектора?

22. Какими командами формируются особые матрицы?

23. Для чего используются символы двоеточие (:),точка с запятой (;) и многоточие (…)?

24. Для чего используются [], ()?

25. Какие функции используются для обработки векторов и матриц?

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

27. Какие команды используются для построения графиков функции двух переменных? В чем их различия?

28. Каков порядок действий для построения графика функции вида y = f(x)?

29. Каков порядок действий для построения графика функции вида z = f(x,y)?

30. Как построить несколько графиков в одних координатных осях?

31. Как можно управлять внешним видом графика?

32. Каким образом можно вывести несколько графиков в разных координатных осях в одном окне?

33. Какие средства предоставляет система ML для построения диаграмм?

34. Как задается полином?

35. Чему равно число элементов в векторе, определяющем полином?

36. Что такое m-файлы?

37. Какие виды m-файлов вы знаете? Чем они отличаются?

38. Каковы правила записи команд в m-файлах?

39. Как создать m-файл?

40. Как вызывается файл-программа в ML?

41. Какова структура файла-функции в ML?

42. Как вызывается файл-функция в ML?

43. Как описать функцию с несколькими входными и выходными параметрами?

44. Какой командой можно ввести данные с клавиатуры в ML?

45. Какой командой можно вывести данные на монитор в ML?

46. Какие операторы цикла существуют в ML?

47. Какие операторы используются для организации ветвлений в ML?

48. Для чего используются функции menu и eval?

49. Для чего используются команды break, return и exit?

50. Какие команды используются для организации диалога в ML?

51. Что означает знак %?