Численное решение линейных дифференциальных уравнений
Рассмотрим уравнение вида
dy(x)/dx = g(x), (11)
где y(x) - неизвестная функция; g(x) - заданная функция.
Уравнение вида (11) называют дифференциальным уравнением первого порядка, поскольку в него входит только первая производная неизвестной функции y = y(x).
В общем случае аналитического решения уравнения (11) не существует. Кроме того, даже если аналитическое решение существует, часто бывает необходимо представить его в графическом виде, чтобы понять его характер.
Эти причины побуждают искать не точные, а приближенные численные решения дифференциальных уравнений.
Типичный метод численного решения дифференциальных уравнений включает в себя преобразование дифференциального уравнения в алгебраическое разностное. Пусть нам необходимо найти решение уравнения (11) в точке х = xk. Положим, что при x = x0 функция y(x) принимает значение y0. Разобьем интервал [x0, xk] на n интервалов шириной Ах. Поскольку уравнение (11) описывает скорость изменения функции y(x) в точке x0, то можно найти приближенное значение функции у в близлежащей точке xi = x0 + Ax, если Ax мало. Будем считать в первом приближении, что функция g(x) постоянна на отрезке [x0, x1]. В этом случае приближенное значение y(x1) определяется выражением
y(xi) = yi ~y(xo) + Ay = y(xo) + g(xo)Ax.
Определив y1, мы можем повторить эту процедуру и найти значение y в точке x2 = x1+Ax:
y(x2) = y2 ~y(xO + Ay = y(x1) + g(x1)Ax.
Очевидным образом это правило можно обобщить и вычислить значение функции в любой точке xt = xi-1+Ax, в том числе в интересующей нас точке xn = xk,, по итерационной формуле
y(xi) = yi ~y(xt-1) + Ay = y(xM) + g(x-1yAx, i = 1, 2, .., n. (12)
Итак, выполнив n шагов вычислений по формуле (12), мы получаем решение уравнения (11).
Задание
Используя метод Эйлера, решите численно дифференциальное уравнение dy/dx = 2x в точке x = 2 c начальными условиями x0 = 1, y0 = 1. Выберите шаг Ax = 0,1. Вычислив приближенное решение y(x), сравните его с точным решением уравнения и вычислите относительную погрешность.
Численное вычисление одномерных интегралов
Рассмотрим определенный интеграл вида
b
F = { f(x)dx
a (13)
Для большинства подынтегральных функций fx) вычислить аналитически данный интеграл не удается, и поэтому интеграл (13) нужно вычислять численно.
Классический метод численного интегрирования основан на геометрической интерпретации интеграла (13) как площади под графиком функции f(x) в пределах от x = а до x = b (рис. 3).
Делим отрезок [a, b] на n равных отрезков длиной Ax:
Ax = (b-a)ln. (14)
Тогда
x0 = a, xi = x0 + iAx, xn = b. (15)
Простейшей оценкой площади под кривой f(x) служит сумма площадей прямоугольников (рис. 4).
В методе прямоугольников значение функции f(x) вычисляется в начале каждого отрезка (в точке слева) и оценка интеграла дается выражением
n-1
F„ =Z f (xi)Ax
<=° (16)
Рис. 3. Геометрический смысл определенного интеграла
Fn =Z f(xi +Ax/2)AX
(17)
Рис. 4. Оценка площади под кривой методом прямоугольников Модификацией метода прямоугольников является вычисление f(x) в средней точке каждого отрезка:
п—1
i =0
Другим приближением является формула трапеций, в которой интеграл оценивается вычислением площади трапеции со сторонами, равными значениям f(x) в начале и конце отрезка (рис. 5).
Рис. 5. Оценка интеграла методом трапеции Оценка интеграла дается выражением
S трапеции = (fxo) + f(xo+dx))/2.
n—1
Fn = [ /(xo)/2 + X f(xi)+ /(xn) ]A*
i=0 (18)
Очевидно, что точность вычисления интегралов вышеописанными методами определяется шагом сетки Лх: чем меньше его величина, тем с большей точностью мы вычисляем интеграл.
Задание
2
Вычислите любым из описанных выше методов интеграл F = J х<^х
i
и сравните приближенный ответ с точным значением интеграла.
Метод Монте-Карло
Метод Монте-Карло является весьма популярным методом в математическом моделировании физических задач и основывается на проведении серии статистических испытаний. Метод Монте -Карло является самостоятельным методом и не относится к числу методов раздела математики, называемого «Численные методы».
Проиллюстрируем работу метода на примере следующей задачи. Каким образом можно измерить площадь озера, имея только кучу камней? Предположим, что озеро расположено на поле известной площади S (рис. 6).
Если бросать камни произвольным образом в пределы озера, то часть камней попадет в озеро. Тогда площадь озера приблизительно равна площади поля, умноженной на долю камней, попавших в озеро. Если n - количество брошенных камней, m - количество камней, попавших в озеро, то площадь озера будет определяться формулой S^epa = S m/n.
Рассмотрим данный метод применительно к вычислению определенного интеграла (рис. 7).
Рис 7. Оценка интеграла методом Монте-Карло Площадь прямоугольника, ограниченного линиями y = H, у = 0, x = a, x = b, равна S. Площадь под кривой у = fx) на интервале [a, b] равна F.
Будем генерировать n пар случайных чисел xt,yt (i = 1, 2, ...n), таких, что:
a < xt < b, 0 (19)
yt < Ax), (20)
Количество m пар точек (xt, yt), которые удовлетворяют условию (20) представляет собой отношение интеграла F к площади прямоугольника S.
F = J f (x)dx
a
Отсюда оценка интеграла методом Монте-Карло определяется выражением Fn = H (b - a) m/n; m - количество точек, удовлетворяющих условию (20); n - общее количество сгенерированных пар точек (xi, yi), удовлетворяющих условию (19).
Другая разновидность метода Монте-Карло основывается на теореме, гласящей, что оценка определенного интеграла определяется следующей
формулой:
b
(22)
F = J f (x)dx & (b-a),
a
где < f > - среднее значение функции f(x) на интервале [a, b].
Для вычисления < f > будем генерировать случайные числа Xi, удовлетворяющие условию a < Xi < b, и вычислять значения функции f(xi) в этих точках. Тогда оценка одномерного интеграла будет иметь вид
b П
(23)
F = J f (x)dx&Fn= (b-a)= (b-a)1 /n- Sf (xi)
^ i=1
где xi - случайные числа; n - число испытаний. Данный метод называется методом выборочного среднего.
Рассмотрим вопрос точности вычисления интегралов методом Монте - Карло. Понятно, что чем больше случайных точек мы сгенерируем, тем точнее будет оценка интеграла, причем можно показать, что точное значение интеграла F лежит в некотором интервале с центром в точке Fn:
Fe (Fn - Gn, Fn + Gn).
Здесь Gn = g/п,
где g - дисперсия; n - число испытаний.
Дисперсия g определяется по формулам:
g2 =2>- 2; (24)
n
(24а)
(246)
<f > = 1/ пЁ f (xi)»
i=1 n
= 1/ nE (/(x))2
Вычисление многомерных интегралов
Многие физические задачи содержат усреднение по многим переменным. Например, предположим, что нам известна зависимость от координаты и скорости полной энергии системы десяти взаимодействующих частиц. Поскольку в трехмерном пространстве каждая частица имеет по три компонента координаты и три компонента скорости, то полная энергия данной системы является функцией шестидесяти переменных. Следовательно, для вычисления средней энергии требуется вычислить шестидесятимерный интеграл. Если разделить область изменения каждой переменной на n отрезков, то в данном случае придется вычислять сумму по n60 точкам. Совершенно очевидно, что для больших размерностей интегралов применять классические численные методы нельзя, тем не менее они используются, когда кратность интегралов не превышает пяти.
Простейший метод оценки многомерных интегралов заключается в сведении этих интегралов к произведению одномерных интегралов.
* * х2 У2
Достарыңызбен бөлісу: |