книги / Несущая способность конструкций в условиях теплосмен

Для всего отрезка [a,b] погрешность интегрирования имеет четвертый порядок относительно шага интегрирования h:

Ψ ≤ M 4 h4 (b − a )720, или Ψ = O(h4 ) .

На рис. 11.20 отображена сходимость приближенного значения определенного интеграла, полученного с помощью формул метода Симпсона.

11.5. Обыкновенные дифференциальные уравнения

Решением, интегралом или интегральной кривой обыкновен-

ного дифференциального уравнения первого порядка

y′ = f (x, y )

x

называется дифференцируемая функция y(x), удовлетворяющая этому уравнению, такая, что y′x (x) ≡ f (x, y (x)) на некотором интерва-

ле изменения аргумента x.

Различают три основных типа дифференциальных задач:

– задача Коши, когда для одного из значений аргумента (например, x0), принимаемого за начальное, известно значение функции y0, то есть y (x0 ) = y0 ;

– граничная задача69, когда условия на искомую функцию задаются для нескольких различных значений аргумента,

y (x0 ) = y0 , y (x1 ) = y1 , y (x2 ) = y2 , ;

– задача на собственные значения, когда в формулировку за-

дачи входят неопределенные параметры, определяемые в ходе решения этой задачи (нахождение частот собственных колебаний многомассовых систем и др.).

В частности, задача Коши заключается в отыскании функции y(x), удовлетворяющей дифференциальному уравнению

69 Граничная задача может быть сформулирована лишь для дифференциальных уравнений второго и более высокого порядков.

171

Здесь f(x, y) – заданная непрерывная функция двух аргументов. Условия существования и единственности решения задачи (11.9)

устанавливает теорема Пеано70 [30]: если функция f(x, y) непрерывна в открытой области D, тогда через каждую точку (u, v) этой области проходит хотя бы одна интегральная кривая. Каждая ин-

тегральная кривая может быть продолжена в обе стороны вплоть до границы любой замкнутой области G, целиком содержащейся в D и содержащей точку (u, v) внутри себя. Кроме того, если функция f(x, y) имеет в области D непрерывные частные производные первого порядка или удовлетворяет условию Липшица

f (x, y2 ) − f (x, y1 ) ≤ K y2 − y1 , (x, y1 ), (x, y2 ) G, K = const > 0,

то задача (10.1) имеет единственное решение.

В дальнейшем предполагается, что правая часть f(x, y) уравнения (11.9) дифференцируема требуемое число раз.

Благодаря развитию современных средств вычислительной техники широкое распространение получили методы численного решения дифференциальных уравнений. Задачи теплопроводности, механики твердого деформируемого тела, жидкостей и газов и многие другие решены благодаря широкому использованию сеточных методов. Ввиду этого понятен интерес к теории разностных схем, позволяющей еще на стадии разработки алгоритма решения сложной инженерной проблемы выяснить условия успешной реализации вычислительной модели.

Численные методы следует применять лишь к корректно поставленным задачам, то есть таким, для которых решение существует, единственно и непрерывно зависит от входных данных, то есть от начальных (граничных) условий и коэффициентов уравнения.

Численное решение задачи считается корректным, если исходная задача поставлена корректно и ее дискретный аналог сохра-

70 Пеано Джузеппе (27.08.1858–20.09.1930) – итальянский математик, профессор Туринского университета.

172

няет свойства корректности, то есть получаемая система алгебраических уравнений однозначно разрешима и устойчива71 по входным данным.

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

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

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

11.5.1. Метод Эйлера

Предполагается, что для интервала [a, b], на котором ищется решение дифференциального уравнения (11.9), построена сеточная область Ω m = {a = x0 < x1 < < xm = b} с постоянным шагом h. Для

построения решения обыкновенного дифференциального уравне-

71 Под устойчивостью в данном случае понимается непрерывная зависимость решения от входных данных [2].

173

ния (11.9) используется разложение искомой функции y(x) в ряд Тейлора вблизи произвольной точки xk Ωm :

y (xk +1 ) = y (xk ) + y′x (xk )h + .

С помощью полученного выражения строится вычислительный процесс

Здесь и далее символами yk обозначается результат численного решения уравнения (11.9), а выражение y (xk ) в дальнейшем используется для обозначения точного решения исходной задачи.

Пример. Решить методом Эйлера дифференциальное уравне-

ние y′ = − y с начальным условием y (0) = 1 .

x

В этом частном случае правая часть дифференциального уравнения имеет вид f (x, y ) = − y , то есть вычислительная схема метода

Эйлера (11.10) записывается следующим образом:

yk +1 = yk − yk h = yk (1 − h), k = 1, 2, .

Строится последовательность значений функции yk для узлов сетки Ωm :

y0 = y (0) = 1;

y1 = y0 (1 − h) = (1 − h);

y2 = y1 (1 − h) = (1 − h)(1 − h) = (1 − h)2 ;

y3 = y2 (1 − h) = (1 − h)3 ; …

ym = (1 − h)m .

174

Таблица 11.10

Результаты численного решения ym методом Эйлера обыкновенного дифференциального уравнения y′ = − y

с начальным условием y(0) = 1

Для оценки погрешности получаемого численного решения производятся вычисления с различными шагами интегрирования h (табл. 11.10). Точное решение рассматриваемой задачи имеет вид y (x) = e− x и, например, для x = 10 принимает значение

y(10) = 0,453 999·10–4.

Три верные значащие цифры получены для шага h = 0,0001, при этом требуется выполнить 100 000 «шагов» по аргументу x, то есть необходимо 100 000 раз выполнить вычисления в соответствии с формулой (11.10).

Погрешность численного решения представляется разностью

zk = yk − y (xk ),

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

Преобразование последней формулы добавлением в правую часть и одновременным вычитанием слагаемого f (xk , y (xk )) приводит к результату

175

+ f (xk , zk + y(xk ))− f (xk , y(xk )).

Выражение в фигурных скобках представляет собой разностный аналог выражения (11.9), в котором численные значения yk искомой функции заменены точными значениями y(xk). Это позволяет проверить, насколько точно разностный аналог (11.11) аппроксимирует исходное дифференциальное уравнение (11.9).

Величина

определяет погрешность аппроксимации разностным аналогом исходного дифференциального уравнения. Для оценки этой погрешности разложение решения в ряд Тейлора

y (xk +1 ) = y (xk ) + y′x (xk )h + y′′xx (xk )h2 2 +

подставляется в формулу для погрешности:

= y′x (xk ) + y′′xx (xk )h2 + − f (xk , y (xk )).

Поскольку, согласно исходному уравнению, y′x (xk ) = f (xk , y (xk )),

оценка погрешности принимает вид

ψ k = y′′xx (xk )h2 + = O (h).

При условии, что значения второй производной y′′xx (xk ) огра-

ничены на интервале [a, b], погрешность аппроксимации дифференциального уравнения (11.9) разностным аналогом (11.11) оказывает-

176

ся величиной, пропорциональной первому порядку шага интегрирования h. В общем случае порядок погрешности аппроксимации равен p, если имеет место соотношение вида

ψk = O(hp ).

На рис. 11.21 представлены результаты решения задачи из рассмотренного примера с использованием формулы (11.10) при малых шагах интегрирования, а также точное решение той же задачи. Видно, что с уменьшением шага интегрирования h численное решение приближается к точному.

0,75

0,5

0,25

0 1 2 3 4 5

в

Рис. 11.21. Сравнение точного решения с результатами расчетов на ЭВМ при различных шагах интегрирования: h = 0,5 (а), h = 0,25 (б), h = 0,1 (в);

— – точное решение; —о— – численное решение

Противоположную картину можно наблюдать при увеличении шага интегрирования. На рис. 11.22 приведены результаты расчетов той же задачи с большими шагами интегрирования.

177

Рис. 11.22. Результаты численного решения разностного уравнения с большими шагами интегрирования: h = 1 (а), h = 1,5 (б), h = 2 (в), h = 2,1 (г); — – точное решение; —о— – численное решение

При шаге h = 1 величины yk уже на втором шаге интегрирования обращаются в 0 и в дальнейшем не изменяются (см. рис. 11.22, а). Шаг интегрирования h = 1,5 дает отрицательные значения yk для некоторых значений xk, что противоречит точному решению задачи (см. рис. 11.22, б). Следует отметить, что даже в этом случае можно наблюдать сходимость численного решения к точному при большом числе шагов. Рис. 11.22, в показывает, что при h = 2 численное решение дает значения yk, осциллирующие возле действительного значения. При дальнейшем увеличении шага получается расходящееся решение (см. рис. 11.22, г).

178

Приведенные результаты говорят о необходимости проведения анализа результатов расчетов, проводимых с использованием ЭВМ, с целью проверки сходимости численных решений к точным.

11.5.2. Метод Рунге – Кутты72 2-го порядка

Для построения разностной схемы используется разложение решения обыкновенного дифференциального уравнения (11.9) в ряд Тейлора:

y (xk +1 ) = y (xk ) + y′x (xk )h + y′′xx (xk )h2 2 + .

Вторая производная вэтом разложении заменяетсявыражением

достижения наибольшей точности записанного выражения.

Для дальнейших выкладок производится замена величины y разложением в ряд Тейлора:

72 Кутта Мартин Вильгельм (03.11.1867–25.12.1944) – немецкий физик и математик, профессор Высшей технической школы Штудгарта.

179

Вводятся обозначения

α = h2 x , β = 1 − h2 x , γ = xh , δ = f (xk , yk ) xh ,

позволяющие записать предыдущее выражение в форме

Очевидно, что введенные коэффициенты зависят от x и могут быть определены через коэффициент α, который в этом случае играет роль параметра:

β = 1 − α, γ = 12α , δ = f (xk , yk )2α .

Окончательно схемаРунге – Куттыпринимаетформу выражения

Та же схема в форме разностного аналога уравнения (11.9) имеет вид

( yk +1 − yk )h = (1 − α ) f (xk , yk ) + αf (xk + h2α , yk + hf (xk , yk )2α ).

При α = 0 получается, как частный случай, известная схема Эйлера

yk +1 = yk + hf (xk , yk ).

При α = 1 выражение (11.12) записывается в форме yk +1 = yk + hf (xk + h2, yk + hf (xk , yk )2).

В этом случае проведение расчетов на очередном шаге интегрирования можно рассматривать как последовательность следующих операций (рис. 11.23):

1. Рассчитывается значение