Рассмотрим расчет переходного процесса в последовательном электрическом колебательном контуре.
Свободные колебания тока в контуре описываются уравнением
L (d²I / dt²) + R (dI / dt) + I / C = 0,
где I - ток, L - индуктивность контурной катушки, R - сопротивление потерь контура, C - емкость конденсатора. Примем R = 5 Ом, L = 10-5 Гн и С = 10-9 Ф.
Используя обозначения z0 = I и z1 = dz0/dt = dI/dt, преобразуем исходное уравнение в систему дифференциальных уравнений первого порядка:
dz0 / dt = z1,
dz1 / dt = -(R z1 + z0/ C) / L,
Приняв в качестве начальных условий I = 1 мА и dI/dt = 0, запишем решение в пакете MathCAD следующим образом:
В этом примере использована функция rkfixed. Она вычисляет решение системы, содержащей обыкновенные дифференциальные уравнения первого порядка методом Рунге-Кутта с фиксированным шагом интегрирования. Обращение к функции выполняется следующим образом:
F:= rkfixed(Z0, t0, tk, N, f),
где F - выходные данные, возвращаемые функцией, Z0 - вектор начальных значений, размерность которого должна соответствовать порядку решаемой системы, t0 - начальное значение аргумента, tk - верхний предел изменения аргумента, N - количество рассчитываемых точек на интервале от t0 до tk, f - имя функции, описывающей правую часть системы дифференциальных уравнений.
Из приведенного примера видно, как записывается функция f(t,z) в документе MathCAD. Имя зависимой переменно z должно совпадать с именами ее компонент, присутствующих в правой части формулы f(t,z). Сами дифференциальные уравнения в предварительно приводятся к каноническому виду (см. систему уравнений выше) так, чтобы в левой части каждого из них находилась только производная по одной из компонент z. Первым в системе должно записываться уравнение, определенное для производной от нулевой компоненты z, то есть от z0, затем - от z1 и так далее. Следует отметить, что отсчет индексов в MathCAD определяется параметром ORIGIN и обычно по умолчанию начинается с нуля.
Функция rkfixed возвращает найденное решение в виде двухмерного массива (см. рис. 1, а). В приведенном примере результат передается переменной F1. Значения независимой переменной t размещаются в нулевом столбце F<0> массива. Они изменяются в пределах от t0 до tk с шагом (tk - t0) / N. В приведенном примере t0 = 0, tk = 5*10-6 c, N = N1 = 500, поэтому шаг изменения аргумента равен (tk -t0) / N1 = 5*10-6 / 500 = 10-8 с.
Найденное численно решение z0(t) и первая производная z1(t) записываются соответственно в первом и во втором столбцах - F1<1> и F1<2>.
Результат решения методом Рунге-Кутта показан графически на рис. 1, б.
Рис. 1. Решение задачи Коши в MathCAD
Выполним аналогичный расчет, уменьшив вдвое шаг при N = N2 = 1000, и оценим погрешность решения по правилу Рунге:
Таким образом, погрешность решения Δ при шаге h/2 не превышает 10-9.
Рассмотрим решение этой же задачи методом Эйлера с использованием рекурсивной процедуры. Изменив для удобства записи обозначения переменных z0 → Y0 и z1 → Y1, преобразуем систему уравнений:
Здесь использована замена производных на их приближенные аналоги dY/dt ≈ ΔY/Δt. Благодаря этому решение задачи можно представить в MathCAD следующим образом:
В данной записи, соответствующей методу Эйлера, параметр Δt - шаг интегрирования, i - счетчик точек решения. Начальные значения переменных t, Y0 и Y1 заданы одновременным присваиванием. Использованная в примере рекурсивная ("возвратная") процедура вычисляет последующие точки решения с индексом i+1 из предыдущих с индексом i.
Рассчитанные методом Эйлера зависимости Y0(t) приведены ниже на рис. 2.
Рис. 2. Расчет методом Эйлера
Зависимость, представленная на рис. 2, а получена при Δt = h = 10-10 с, а на рис. 5, б - при Δt = h = 7 * 10-9 с. В последнем случае неприемлемый с физической точки зрения результат объясняется выбором слишком большого шага интегрирования (см. п. 1.2.4). Метод Эйлера дает корректный результат при h = 10-10 с, то есть расчете порядка 50000 точек. Однако оценка погрешности такого решения по правилу Рунге дает Δ ≈ 1,5 * 10-3.
Кроме rkfixed в MathCAD имеется еще ряд функций, обеспечивающих решение задачи Коши. К ним относятся odesolve, rkadapt, Bulstoer, Stiffb и Stiffr. Функциии rkadapt и Bulstoer предназначены для решения задач с гладкими функциями и медленно меняющимися решениями, Stiffb и Stiffr - для задач с жесткими системами. Функция sbval позволяет привести краевую задачу к задаче Коши и решить ее методом стрельбы.