Численное решение дифференциальных уравнений

Лекция No 6 Численное решение дифференциальных уравнений Рассмотрим стандартные средства численного решения дифференциальных уравнений и систем дифференциальных уравнений. 1. Численное решение одного дифференциального уравнения MathCAD 2000 позволяет без дополнительных преобразований численно решить дифференциальное уравнение, явно разрешенное относительно старшей производной (рис. 15). а) б) Рис. 15. Примеры уравнений, разрешенного (а) и неразрешенного (б) относительно старшей производной Решение осуществляется с помощью специального блока Given-Odesolve, состоящего из следующих компонент: 1. Директива Given. 2. Дифференциальное уравнение, записанное в традиционной математической форме со следующими особенностями: а) вместо простого знака равенства «=» используется оператор логического равенства (вводится нажатием Ctrl-=); б) при обозначении интегрируемой функции всегда указывается аргумент (то есть вместо функции x(t) нельзя писать просто x); в) при записи производных используются либо стандартные операторы и , либо ставятся (с помощью Ctrl-F7) символы производной, например x’(t), x’’(t). 3. Указание начальных или конечных значений интегрируемой функции и ее производных (за исключением старшей), входящих в уравнение. Значения вводятся в традиционной форме с использованием оператора логического равенства. Число значений должно совпадать с порядком уравнения. Для уравнения второго порядка вида должны быть заданы начальные значения функции и ее первой производной, например x(0) = 1; x’(0) = 0,5. Для ввода символа производной «’» используется комбинация клавиш Ctrl-F7. 4. Обращение к функции Odesolve. Первый аргумент - всегда имя независимой переменной. Второй аргумент - конечное значение независимой переменной. Третий (необязательный) аргумент - количество промежуточных точек решения. Odesolve возвращает функцию, представляющую приближенное (численное) решение дифференциального уравнения на заданном интервале времени. Данная функция может быть использована для определения значений интегрируемой функции в различных точках, а также для построения графика. Пример. Решим вышеуказанное дифференциальное уравнение при значениях t = 0..5; найдем значения x при t = 2; 4, и построим график решения. Численное решение систем дифференциальных уравнений с использованием Odesolve возможно только начиная с версии MathCAD 11. В более старых версиях можно воспользоваться специальными функциями (rkadapt, rkfixed, bulstoer). 2. Численное решение систем обыкновенных дифференциальных уравнений Дифференциальные уравнения, входящие в систему, должны иметь первый порядок (то есть содержать только первые производные). Все уравнения должны быть предварительно разрешены относительно производных и записаны в нормальной форме вида . Для преобразования уравнений в нормальную форму есть два основных подхода: 1. Понижение порядка уравнений путем замены переменных. Если исходное дифференциальное уравнение порядка q (q>1) имеет вид , то вводятся новые переменные pj, причем j = 1..q-1. В исходном уравнении производится серия замен: , а производная высшего порядка заменяется производной первого порядка: . Добавляется q – 2 новых уравнений вида . Добавляется еще одно уравнение . Например, уравнение можно преобразовать в систему уравнений: 2. Приведение системы дифференциальных уравнений к явному виду путем ее решения относительно производных. Например, решая систему относительно и , получим: Рассмотрим решение систем дифференциальных уравнений в MathCAD на примере задачи о моделировании динамики электрической цепи, показанной на рис. 16. Динамика описывается следующей системой дифференциальных уравнений: Рис. 16. Электрическая цепь где Uc - напряжение на конденсаторе. Пусть ; i1(0) = i2(0) = Uc(0) = 0; ; L1 = 0,02; L2 = 0,06; M = 0,01; R = 0,5; C = 0,01. Решение записывается следующим образом: 1. Определяются все константы и вспомогательные функции, присутствующие в правой части системы. 2. Определяется специальная функция, вычисляющая правую часть системы. Функция имеет два аргумента: первый - независимая переменная (например, время t), второй - вектор текущих значений зависимых переменных. Результатом функции должен быть вектор, содержащий значения правых частей системы, вычисленных по значениям второго аргумента функции. Векторы имеют столько элементов, сколько уравнений в системе. При записи правых частей все зависимые переменные заменяются элементами вектора – второго аргумента, причем используется следующее правило: нулевому элементу соответствует переменная, производная от которой стоит в левой части первого уравнения; первому элементу - переменная, производная от которой стоит в левой части второго уравнения и так далее. В приведенном далее примере, где второй аргумент функции обозначен как Y, элементу Y0 соответствует i1 - переменная из производной в левой части первого уравнения, элементу Y1 соответствует i2 - переменная из производной в левой части второго уравнения, элементу Y2 соответствует UC. 3. Задается вектор начальных значений независимых переменных. 4. Обращение к функции rkfixed. Первый аргумент - вектор начальных значений. Второй и третий - соответственно начальное и конечное значения независимой переменной. Четвертый аргумент - число промежуточных точек решения (обычно достаточно большое число в диапазоне ). Пятый - имя функции, вычисляющей правую часть системы. Функция rkfixed возвращает матрицу, в нулевом столбце которой находятся значения независимой переменной, а в прочих столбцах соответствующие значения зависимых переменных. Решение показано на рис. 17. Рис. 17. Запись решения задачи в MathCAD На рис. 18 показаны графики i2(t), Uc(t). Данным переменным соответствуют второй и третий столбцы матрицы S. Рис. 18. Графики i2(t), Uc(t)