Численное решение одного дифференциального уравнения

Численное решение одного дифференциального уравнения

а)                б)

Рис. 15. Примеры уравнений, разрешенного (а) и неразрешенного (б) относительно старшей производной

Решение осуществляется с помощью специального блока Given-Odesolve, состоящего из следующих компонент:

2. Дифференциальное уравнение, записанное в традиционной математической форме со следующими особенностями: а) вместо простого знака равенства «=» используется оператор логического равенства (вводится нажатием Ctrl-=); б) при обозначении интегрируемой функции всегда указывается аргумент (то есть вместо функции x(t) нельзя писать просто x); в) при записи производных используются либо стандартные операторы и , либо ставятся (с помощью Ctrl-F7) символы производной, например x’(t), x’’(t).

3. Указание начальных или конечных значений интегрируемой функции и ее производных (за исключением старшей), входящих в уравнение. Значения вводятся в традиционной форме с использованием оператора логического равенства. Число значений должно совпадать с порядком уравнения. Для уравнения второго порядка вида должны быть заданы начальные значения функции и ее первой производной, например x(0) = 1; x’(0) = 0,5. Для ввода символа производной «’» используется комбинация клавиш Ctrl-F7.

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

Численное решение систем дифференциальных уравнений с использованием Odesolve возможно только начиная с версии MathCAD 11. В более старых версиях можно воспользоваться специальными функциями (rkadapt, rkfixed, bulstoer).

2. Численное решение систем обыкновенных дифференциальных уравнений

Для преобразования уравнений в нормальную форму есть два основных подхода:

1. Понижение порядка уравнений путем замены переменных. Если исходное дифференциальное уравнение порядка q (q>1) имеет вид

,

то вводятся новые переменные p_j, причем j = 1..q-1. В исходном уравнении производится серия замен:

а производная высшего порядка заменяется производной первого порядка:

Добавляется q – 2 новых уравнений вида

Добавляется еще одно уравнение

Например, уравнение

можно преобразовать в систему уравнений:

2. Приведение системы дифференциальных уравнений к явному виду путем ее решения относительно производных. Например, решая систему

относительно и , получим:

Рассмотрим решение систем дифференциальных уравнений в MathCAD на примере задачи о моделировании динамики электрической цепи, показанной на рис. 16. Динамика описывается следующей системой дифференциальных уравнений:

Рис. 16. Электрическая цепь

где U_c - напряжение на конденсаторе. Пусть ; i₁(0) = i₂(0) = U_c(0) = 0; ; L₁ = 0,02; L₂ = 0,06; M = 0,01; R = 0,5; C = 0,01.

1. Определяются все константы и вспомогательные функции, присутствующие в правой части системы.

2. Определяется специальная функция, вычисляющая правую часть системы. Функция имеет два аргумента: первый - независимая переменная (например, время t), второй - вектор текущих значений зависимых переменных. Результатом функции должен быть вектор, содержащий значения правых частей системы, вычисленных по значениям второго аргумента функции. Векторы имеют столько элементов, сколько уравнений в системе. При записи правых частей все зависимые переменные заменяются элементами вектора – второго аргумента, причем используется следующее правило: нулевому элементу соответствует переменная, производная от которой стоит в левой части первого уравнения; первому элементу - переменная, производная от которой стоит в левой части второго уравнения и так далее. В приведенном далее примере, где второй аргумент функции обозначен как Y, элементу Y₀ соответствует i₁ - переменная из производной в левой части первого уравнения, элементу Y₁ соответствует i₂ - переменная из производной в левой части второго уравнения, элементу Y₂ соответствует U_C.

3. Задается вектор начальных значений независимых переменных.

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

Решение показано на рис. 17.

Рис. 17. Запись решения задачи в MathCAD

На рис. 18 показаны графики i₂(t), U_c(t). Данным переменным соответствуют второй и третий столбцы матрицы S.

Рис. 18. Графики i₂(t), U_c(t)

 

Некоторые стандартные функции MathCAD

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

M - квадратная матрица;

V - вектор (матрица из одного столбца);

A - произвольная матрица;

S - симметричная матрица;

G - произвольная матрица или число;

X - вектор или число;

Z - комплексная матрица или число;

z - комплексное число;

прочие символы - скалярные величины.

 

Экспоненциальные и логарифмические функции

exp(X) - экспонента от X;

ln(X) - натуральный логарифм от X;

log(X) - десятичный логарифм от X;

log(X,b) - логарифм от X по основанию b.

 

Гиперболические и тригонометрические

sin(X), cos(X), tan(X), cot(X), sec(X), csc(X) - соответственно синус, косинус, тангенс, котангенс, секанс, косеканс от X, причем аргументы указываются в радианах;

sinh(X), cosh(X), tanh(X), coth(X), sech(X), csch(X) - аналогичные гиперболические функции;

asin(z), acos(z), atan(z), acot(z), asec(z), acsc(z) - соответственно арксинус, арккосинус, арктангенс, арккотангенс, арксеканс, арккосеканс от z.

 

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

Re(Z), Im(Z) - соответственно вещественная и мнимая части комплексного числа Z;

arg(z) - аргумент комплексного числа z (в радианах).

 

Матричные функции

length(V) - возвращает число элементов вектора V;

cols(A) - возвращает число столбцов матрицы A;

rows(A) - возвращает число строк матрицы A;

matrix(m,n,f) - матрица размером mxn, значения элементов матрицы определяются f - функцией f(i,j) от двух переменных (номера строки и номера столбца). Эта функция должна быть предварительно определена пользователем;

identity(n) - единичная матрица ;

tr(M) - след матрицы M (сумма элементов главной диагонали);

rank(A) - ранг матрицы M;

norme(M) - эвклидова норма матрицы M, то есть корень квадратный из суммы квадратов всех элементов;

eigenvals(M) - вектор, элементы которого являются собственными числами матрицы M;

eigenvecs(M) - матрица, состоящая из нормализованных собственных векторов матрицы M;

cholesky(S) - возвращает нижнетреугольную матрицу L - результат разложения Холецкого вида ;

lu(M) - возвращает матрицу размера , состоящую из трех соединенных матриц P, L, U, являющихся результатом LU-разложения вида .

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

 

Элементы статистического анализа данных

gmean(G1,G2,G3…) - среднее геометрическое аргументов;

mean(G1,G2,G3…) - среднее арифметическое аргументов;

var(G1,G2,G3…) - дисперсия;

stdev(G1,G2,G3…) - среднеквадратичное отклонение.

 

Дискретные преобразования

fft(V1), ifft(V2) - прямое и обратное быстрые преобразования Фурье над вещественными данными. V1 - вектор из 2^m элементов, V2 - вектор из 1 + 2^m-1 элементов, m>2;

cfft(A), icfft(A) - прямое и обратное преобразования Фурье над вещественными и комплексными векторами и матрицами;

wave(V), iwave(V) - прямое и обратное вейвлет-преобразования, V - вектор из 2^m элементов, m - целое число.

 

Аппроксимация, интерполяция и экстраполяция

Аппроксимация - поиск функции, которая с заданной степенью точности описывает исходные данные.

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

Экстраполяция - определение наиболее правдоподобных последующих значений на основании анализа предыдущих значений (предсказание дальнейшего поведения неизвестной функции).

Применяются следующие функции MathCAD:

regress(VX,VY,k) - возвращает вектор данных, используемый для поиска интерполирующего полинома порядка k. Полином должен описывать данные, состоящие из упорядоченных значений аргумента (VX) и соответствующих значений неизвестной функции (VY), то есть график полинома должен проходить через все точки, заданные координатами (VX,VY), или максимально близко к этим точкам;

interp(VS,VX,VY,x) - возвращает интерполированное значение неизвестной функции при значении аргумента x. VS - вектор значений, который вернула функция regress. VX,VY - те же данные, что и для regress. Функции interp и regress используются в паре;

predict(V,m,n) - возвращает вектор из n предсказанных значений на основании анализа m предыдущих значений из вектора V. Предполагается, что значения функции в векторе V были получены при значениях аргумента, взятых последовательно, с одинаковым шагом. Используется алгоритм линейной предикции. Наиболее целесообразно использовать predict для предсказания значений по данным, в которых отмечены колебания.

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

1. Пусть заданы координаты пяти точек (1; 1), (2; 2), (3; 3), (4; 2), (5; 3), представляющих результаты измерения значений некоторой неизвестной функции при различных значениях x. Необходимо подобрать интерполирующую функцию (гладкую кривую), проходящую через заданные точки.

2. Дана функция y(i) = e^–i/10*sin(i). Известны значения данной функции при i = 0, 1, …, 10. Основываясь на десяти последних значениях, необходимо предсказать последующие десять значений.

Решения показаны на рис. 19.

а)                    б)

Рис. 19. Решения в MathCAD первой (а) и второй (б) задач

 

Нахождение корней полинома

polyroots(V) - возвращает вектор, содержащий все корни полинома , заданного вектором-столбцом коэффициентов

.

 

Прочие функции

max(G1,G2,…) - максимальное значение среди аргументов;

min(G1,G2,…) - минимальное значение среди аргументов;

if(a,b,c) - возвращает b, если , иначе возвращает c;

sign(a) - возвращает –1, 0 или 1 в зависимости от знака числа a.

На рис. 20 показан пример применения функции if.

Рис. 20. Функция, вычисляющая факториал

 

Элементы программирования в MathCAD

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

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

Программный блок представляет собой группу операторов присваивания и управляющих операторов. Необходимо обратить особое внимание, что все ключевые слова (например, if) в этих операторах обязательно вводятся с помощью панели Programming (Программирование). Их ввод с клавиатуры - ошибка!

В целом правила работы с операторами те же, что и в языке Pascal, отличия касаются способа записи операторов.

Таблица 2. Соответствие программных операторов MathCAD и Pascal

Оператор языка Pascal	Оператор MathCAD	Комментарий
A := B		Присваивание
Begin оператор1; оператор2; … End		Группа, объединяющая несколько операторов в один составной оператор. Для создания группы и добавления в нее новой пустой строчки используется кнопка «Add Line» панели Programming
If условие Then оператор   If условие Then Begin оператор1; оператор2; … End	оператор if условие	Простой оператор ветвления. Как и в языке Pascal, его действие распространяется на один указанный оператор, который может быть группой операторов. Условием может быть любое логическое выражение, которое может содержать знаки отношения (вместо обычного знака равенства используется знак логического равенства) и логические операторы (находятся на панели Boolean): - Not; - And; - Or; - Xor
If условие Then оператор1 Else оператор2		Полный оператор ветвления
For инд := нач To кон Do оператор		Фиксированный оператор цикла. Индексная переменная принимает значения от начального до конечного с шагом, равным единице. Цикл действует на один указанный оператор, который может быть группой операторов
While условие Do оператор		Гибкий оператор цикла с предусловием. Цикл выполняется, пока истинно заданое условие
Break Continue	break continue	Оператор break принудительно завершает текущий цикл. Оператор continue завершает только текущий виток цикла и начинает следующий виток
Нет прямого аналога	выражение1 on error выражение2	Специальная операция обработки ошибок. Сначала вычисляется выражение2. Если при этом происходит ошибка, то результатом операции будет выражение1. Если ошибки нет, то результат - выражение2. Пример: Здесь локальная переменная A получает значение 2, переменная B - значение 0,5

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

return значение,

который также вводится с помощью панели Programming. Его выполнение заканчивает работу текущего программного блока.

1. Функция, возвращающая –1, 0 или 1 в зависимости от знака аргумента.

2. Пусть интегрируется дифференциальное уравнение

где параметр z определяется в результате решения нелинейного уравнения

Известно, что в рассматриваемом случае это нелинейное уравнение имеет единственное решение. Создадим функцию, которая решает данное уравнение методом касательных с заданной точностью ?.

Функция Solve возвращает значение z, которое является корнем уравнения при заданном значении x. Решение дифференциального уравнения: