3.3.1. Повышение порядка точности аппроксимации обыкновенных дифференциальных уравнений
Однако порядок точности аппроксимации дифференциальных уравнений можно существенно увеличивать, например, вводя полуцелые узлы.
Для простоты рассуждений рассмотрим следующее обыкновенное дифференциальное уравнение с начальным условием:
;
N(t0)
= N0;
.
(3.3.1)
Отрезок [t0, tk] разобьем на k равных частей, при этом выполняется следующая цепочка неравенств:
0 = t0 t1 . . .< tj <… tk - 1 tk = tk, (3.3.2)
где tj = j τ; τ = ( tk – t0 ) / k; j = 0, … , k .
Алгоритм численного решения уравнения (3.3.1) на разностной сетке (3.3.2) с введением полуцелых узлов состоит из двух шагов: на первом полушаге по времени τ / 2 вычисляется численное решение уравнения с первым порядком точности и на втором полушаге уточнения численного решения со вторым порядком точности. Для записи дифференциального уравнения (3.3.1) в разностном виде на этих шагах рассмотрим разностную схему, приведенную на рис. 3.22.
N
Nj
Nj+1
tk
t
0
t0
j
j+1
k
Рис. 3.22. Разностная схема для решения дифференциального уравнения (3.3.1)
На
первом полушаге по времени τ / 2 (на
разностной сетке вводились полуцелые
временные слои
)
аппроксимация дифференциального
уравнения (3.3.1) проводилась на основе
правой разностной производной, имеющий
первый порядок точности, следующим
образом (схема Эйлера):
или
;j
= 0, …, k-1.
(3.3.3)
На втором полушаге по времени τ / 2 (этап – коррекции) необходимо исходное уравнение аппроксимировать со вторым порядком точности, используя центральную разностную производную:
или
;j
= 0, …, k-1.
(3.3.4)
[III, 10-11, 17-19].
Результаты численных расчетов. Для практического определения роли порядка точности аппроксимации дифференциального уравнения, при его численном решение, графически сравнивались результаты расчетов, полученных по разностным схемам, имеющих второй порядок точности (3.3.3 – 3.3.4) и по схемам с первым порядком точности:
или
;j
= 0, …, k-1,
(3.3.5)
при
задании следующих исходных данных:
;
tk
= 3600c;
k
= 3000; N0
= 3.
Текст программы алгоритма решения обыкновенного дифференциального уравнения (3.3.1), написанной на блочно-структурном языке системы MATH CAD и числовые данные приведены ниже.
Результаты численных расчетов функции N(t) представлены на рис.3.23, здесь и далее сплошной кривой отмечено решение уравнения (3.3.1), полученное на основе разностных схем, имеющих второй порядок точности, а пунктирной - первый порядок точности.
Рис.3.23. Результаты численных расчетов функции N(t) при шаге
интегрирования τ = 12с и τ = 120с
Как следует из анализа рис. 3.23, максимальное расхождение значений функции N(t) при временном шаге τ = 12с, полученных на основе двух разностных схем, составляет порядка десяти процентов и, соответственно, при таких шагах τ могут использоваться разностные схемы, имеющие первый порядок точности.
Таким образом, обе разностные схемы счетно-устойчивы, но разностные схемы, имеющие первый порядок точности, при больших значениях временных шагов τ приводят к значительным погрешностям (см. рис.3.23), делающих непригодными эти схемы при интегрировании уравнений. В то же время разностные схемы, имеющие второй порядок точности, являются вполне приемлемыми, так как погрешность численного решения не более 15% (см. рис. 3.23).
Задание. Задать нелинейную функцию правой части уравнения (3.3.1) f(N, t) и определить максимальный шаг интегрирования τ, при котором относительная погрешность вычисления функции N(t), полученная на основе разностных схем, имеющих первой и второй порядки точности, не превышает 5%.