Численное решение систем обыкновенных дифференциальных уравнений
Дифференциальные уравнения, входящие в систему, должны иметь первый порядок (то есть содержать только первые производные). Все уравнения должны быть предварительно разрешены относительно производных и записаны в нормальной форме вида
.
Для преобразования уравнений в нормальную форму есть два основных подхода:
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)
Лекция No 7
Некоторые стандартные функции 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 - вектор из 2m элементов, V2 - вектор из 1 + 2m-1 элементов, m>2;
cfft(A), icfft(A) - прямое и обратное преобразования Фурье над вещественными и комплексными векторами и матрицами;
wave(V), iwave(V) - прямое и обратное вейвлет-преобразования, V - вектор из 2m элементов, 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. Функция, вычисляющая факториал
Лекция No 8
Элементы программирования в 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. Решение дифференциального уравнения:
