Лекция №4 ММвХТ 8 ноября 12
Решение обыкновенных дифференциальных уравнений в smath studio (задача коши) Введение
Инженерные исследования динамики процессов, протекающих в механизмах, реакторах, локальных системах стабилизации параметров технологических процессов, трубопроводах, теплообменных процессах и других химических объектах приводят к дифференциальным уравнениям, т.е. уравнениям, содержащим производные. Например: уравнение колебания массы под воздействием силы, изменение концентрации растворов, уравнение накопления и стекания заряда на обкладках конденсатора.
Классические методы решения обыкновенных дифференциальных уравнений часто на практике либо приводят к сложным решениям, либо вообще неприменимы (коэффициенты или функции в дифференциальном уравнении содержат существенные нелинейности либо заданы в виде таблиц экспериментальных данных).
Поэтому большое значение приобретают методы приближенного интегрирования дифференциальных уравнений, которые в зависимости от формы представления решения можно разделить (условно) на три группы:
1) аналитические, дающие приближенное решение дифференциального уравнения в виде аналитического выражения;
2) графические, дающие приближенное решение в виде графика;
3) численные, дающие приближенное решение в виде таблицы.
Рассмотрим численные методы решения обыкновенных дифференциальных уравнений.
1. Постановка задачи численного метода решения дифференциальных уравнений
Дано обыкновенное дифференциальное уравнение первого порядка
у' = f (х, у). (1)
Тре6уется найти решение y=y(х) этого уравнения, удовлетворяющее начальному условию
у(х0) == y0 . (2)
Такая задача называется задачей Коши. Геометрический смысл ее решения состоит в нахождении интегральной кривой у=у(х), проходящей через заданную точку А0 (х0,у0) (рис. 1).
Существуют следующие постановки задачи Коши:
дифференциальное уравнение первого порядка, разрешённое относительно производной:
y ' = f(x,y)
y(x0) = y0,
система n дифференциальных уравнений первого порядка, разрешённая относительно производных (нормальная система n-го порядка):
дифференциальное уравнение n-го порядка y(n)=f(x,y,y',…y(n-1)), разрешённое относительно старшей производной y(n):
Численное решение задачи Коши состоит в нахождении значений y1, y2,...,уn в точках х1=х0+h, х2=х0+2h,…, хn=хо+nh отрезка [a, b], где h - шаг интегрирования, х0=а, хn=b.
|
|
Рис. 1. График кривой у=у(х) |
Рис. 2. Ломаная Эйлера |
Нанеся точки (хо,уо), (х1,у1),...,(хn,уn) на координатную плоскость и соединив отрезками прямой, получим ломаную линию, называемую ломаной Эйлера, - приближенное изображение интегральной кривой (рис.2).
Самый простой и менее точный - метод Эйлера (метод первого порядка точности, расчетные формулы (4)). Наиболее распространенным является метод Рунге-Кутта (метод четвертого порядка точности, расчетные формулы (5)).
Последовательность вычислений по указанным методам следующая:
Разбиваем отрезок [а,b] на n равных частей точками хi=хо+ih (i = 1. 2. ..., n), где шаг h =(b-a)/n, хо=а, хn=b.
Находим решение уравнения, рассчитывая для каждого i (i=0,1,2,…,n) значения
yi+1 = yi + Δyi, (3)
где Δyi= k1(i) (4)
для метода Эйлера и
Δyi=1/6(k1(i)+2 k2(i)+ 2k3(i) + k4(i)) (5)
для метода Рунге-Кутта.
Коэффициенты kj(i), j=1,4 для выражений (4) и (5) определяются по следующим формулам:
k1(i)=hf(xi,yi);
k2(i)=hf(xi+h/2, yi+ k1(i)/2);
k3(i)=hf(xi+h/2, yi+ k2(i)/2);
k4(i)=hf(xi+h, yi+ k3(i)).
Для выполнения вычислений по методу Рунге-Кутта вручную удобно пользоваться схемой, приведенной в табл. 33. Ориентировочную оценку погрешности метода Рунге-Кутта ε можно вычислить по формуле (6)
ε=|y2h - yh|/15. (6)
Метод Рунге-Кутта является одним из методов повышенной точности и, несмотря на свою трудоемкость, широко используется при численном решении, как дифференциальных уравнений, так и систем дифференциальных уравнений.
Среди встроенных функций SMath Studio для численного решения дифференциальных уравнений есть встроенная функция, реализующая метод Рунге – Кутта с постоянным фиксированным шагом. Она имеет вид: rkfixed( v,x0,xk,n,F). Здесь v - начальные условия, записанные в виде вектора, x0, xk – начальное и конечное значения аргумента, n- число шагов, F- правые части системы, записанные в виде вектора. Длина шага интегрирования определяется выражением h=(xk-x0)/n.
Как и большинство методов численного решения дифференциальных уравнений, метод Рунге – Кутта требует предварительного представления дифференциального уравнения n-го порядка в виде системы дифференциальных уравнений первого порядка:
Преобразование дифференциального уравнения в систему уравнений будет рассмотрен при решении примеров.
Отметим, что в последних версиях MathCad появилась функция odesolve(х, b) (ordinary differential equation solution – решение обыкновенного дифференциального уравнения), позволяющая решать уравнение без его преобразования. Здесь х – переменная интегрирования, b- верхняя граница изменения аргумента. Нижняя граница равна нулю.