Численные методы

Во многих научных и инженерных задачах возникает необходимость решения уравнений вида:

f(x, p₁, p_2, , … p_n) = 0 (1)

где f - заданная функция; x - неизвестная величина;, p₁, p_2, , … p_n - параметры задачи.

Функция f(x) может быть алгебраической или трансцендентной. Будем предполагать, что она дифференцируема.

В общем случае аналитических формул для нахождения корней нет и приходится пользоваться приближенными методами нахождения корней.

Они состоят из двух этапов:

1) отыскание первоначального (грубого) приближенного значения корня;

2) уточнение приближенного значения до заданной степени точности.

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

На втором этапе производится последовательное, шаг за шагом, уточнение первоначального грубого приближения. Каждый такой шаг уточнения называется итерацией, а численный метод - методом итераций.

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

Метод простых итераций

От исходного уравнения f(x)= 0 переходят к эквивалентному уравнению:

x = j(x) (2)

Пусть первоначальное приближение будет равно x=x ₀. Подставляем его в правую часть уравнения (2) и получаем новое приближение x ₁ = j(x ₀ ), затем x ₂ = j(x ₁ ), и так далее:

x _k+1 = j(x _k ), (3)

Итерационный процесс будет сходиться к корню уравнения, если выполняется условие:

│ j^’(x)│ < 1 (4)

Чем меньше │ j^’(x)│ вблизи корня, тем быстрее сходится процесс. Графическое пояснение метода простых итераций приведено на рис.1.

Переход от уравнения (1) к уравнению (2) можно осуществить различными способами в зависимости от вида функции f(x). Рассмотрим один из алгоритмов такого перехода.

Умножим левую и правую части уравнения (1) на произвольную константу m и добавим к обеим частям неизвестное x

x + m × f(x) = x.

Введем обозначение:

j(x) = x + m × f(x) (5)

Произвольный выбор константы m позволяет обеспечить выполнение условия сходимости │ j^’(x)│ < 1.

Желательно выбирать значение m таким, чтобы соблюдалось неравенство -1 < j^’(x) < 0. Тогда сходимость итерационного процесса будет двухсторонней. Критерием окончания итерационного процесса является неравенство │ x_k+1 – x_k│ < e, где e - заданная абсолютная погрешность вычисления корня.

Если функция j(x) выбрана по формуле (5), то имеем:

j^’(x) = 1 + m× f ^’ (x).

Наибольшую скорость сходимости получим при j^’(x) = 0, тогда

m = - 1/ f ^’ (x).

В этом случае имеем итерационную формулу Ньютона-Рафсона:

x _k+1 = x _k - f(x _k )/f ^‘ (x _k ) (6)

При этом метод итераций называют методом касательных. На рис. 2 приведено графическое пояснение метода касательных.

Рассмотрим пример численного решения уравнения. На торце полубесконечного стержня (рис. 3) поддерживается постоянная температура T _m.

Рис.1. Графическое пояснение метода простых итераций:

а – односторонний сходящийся процесс (0 < j^’(x) < 1)

б – односторонний расходящийся процесс (j^’(x) > 1)

в – двухсторонний сходящийся процесс (-1 < j^’(x) < 0)

г – двухсторонний расходящийся процесс (j^’(x) < -1)

Рис. 2. Графическое пояснение метода Ньютона-Рафсона.

Рис. 3. Температурное поле в полубесконечном стержне при

постоянной температуре на его торце: T( 0, t)=T _m

Составить программу для расчета глубины x прогрева стержня до температуры T в заданный момент времени t.

Температурное поле в стержне записывается уравнением: (7)

Здесь a - коэффициент температуропроводности металла стержня.

f(x) = T - T(x, t). Используя формулу Ньютона-Рафсона, имеем:

Ниже приведены программы расчета методом простых итераций и методом касательных.

В строках 10-30 производится расчет функции erf(u). При этом используется ее приближение полиномом:

erf(u) = 1 -(a ₁ t + a ₂ t ² +...+ a ₅ t ⁵ ) × exp(-u), где t = 1 /( 1 +px).

Далее приведены результаты расчета по этим программам глубины прогрева до температуры T =500 в момент времени t =5 при T _m=1000 и a =0.1.

При расчете методом касательных число итераций в 3 раза меньше.

Программа решения уравнения методом простых итераций:

1 PRINT " Решение уравнения методом простых итераций"

5 PI=3.14159265358979

10 DEF FNA(X)=1/(1+.3275911*X)

20 DEF FNB(X)=.254829592*FNA(X)-.284496736*FNA(X)^2+1.421413741*

FNA(X)^3-1.453152027*FNA(X)^4+1.061405429*FNA(X)^5

30 DEF FNE(X)=1-FNB(X)*EXP(-X^2)

40 INPUT " Tm, T, t, a ="; TM, TA, S, A

50 X=0: H=1: M=0

60 X=X+H: M=M+1: IF M=50 GOTO 140

70 T=TM*(1-FNE(X/2/SQR(A*S)))

80 IF T< =TA GOTO 100

90 GOTO 60

100 IF H< 0.00001 GOTO 120

110 X=X-H: H=H/2: GOTO 60

120 PRINT " M=" M; " X=" X

130 GOTO 40

140 PRINT " M=" M; " Процесс итераций остановлен"

150 END

Расчет на ЭВМ:

Решение уравнения методом простых итераций

Tm, T, t, a =? 1000, 500, 5, 0.1

M= 27 X=.6744918823242188

Tm, T, t, a =?

Программа решения уравнения методом итераций Ньютона-Рафсона:

1 PRINT " Решение уравнения методом итераций Ньютона-Рафсона"

5 PI=3.14159265358979

10 DEF FNA(X)=1/(1+.3275911*X)

20 DEF FNB(X)=.254829592*FNA(X)-.284496736*FNA(X)^2+1.421413741*

FNA(X)^3-1.453152027*FNA(X)^4+1.061405429*FNA(X)^5

30 DEF FNE(X)=1-FNB(X)*EXP(-X^2)

40 INPUT " Tm, T, t, a ="; TM, TA, S, A

50 X=1: E=0.00001: M=0

60 GOSUB 90: X=X-F: M=M+1

65 IF M=50 GOTO 100

70 IF ABS(F)> E THEN 60

80 PRINT " M=" M; " X=" X: GOTO 40

90 F=(TA-TM*(1-FNE(X/2/SQR(A*S))))*SQR(PI)/TM/2/EXP(-(X^2/4/A/S)): RETURN

100 PRINT " M=" M; " Процесс итераций остановлен"

120 END

Расчет на ЭВМ:

Решение уравнения методом итераций Ньютона-Рафсона

Tm, T, t, a =? 1000, 500, 5, 0.1

M= 9 X=.6744924783706665

Tm, T, t, a =?

Необходимо решить систему из n линейных алгебраических уравнений (СЛАУ):

a ₁₁ x ₁ + a ₁₂ x ₂ +... + a _1n x _n = b ₁

a ₂₁ x ₁ + a ₂₂ x ₂ +... + a _2n x _n = b ₂ (8)

..................

a _n1 x ₁ + a _n2 x ₂ +... + a _nn x _n = b _n

Часто используется матричная форма записи этой системы:

AX = B,

где A - матрица коэффициентов, X - вектор неизвестных, B - вектор свободных членов.

где i – номер строки (i = 1, 2, …n), j – номер столбца (j = 1, 2, …n).

Для решения СЛАУ используются в основном два класса методов: прямые (точные) и итерационные.

Прямые методы являются универсальными и применяются для систем сравнительно невысокого порядка (n< 200).

Итерационные методы выгодно использовать для СЛАУ высокого порядка со слабо заполненными матрицами коэффициентов.

К прямым методам относится метод Гаусса (метод исключения).

Алгоритм метода Гаусса состоит из двух этапов. Первый этап называется прямым ходом метода и заключается в последовательном исключении неизвестных из уравнений, начиная с x ₁.

Из первого уравнения системы выражаем неизвестное x ₁:

x ₁ = (b ₁ - a ₁₂ x ₂ -... – a _1n x _n )/a ₁₁.

Это выражение подставляем во все остальные уравнения системы. Тем самым исключаем x ₁ из всех уравнений, кроме первого.

Затем неизвестное x ₂ выражаем из второго уравнения системы и исключаем его из всех последующих уравнений и т.д.

В результате получаем треугольную систему уравнений:

Эта система имеет верхнюю треугольную матрицу коэффициентов, у которой все элементы ниже главной диагонали равны нулю.

Второй этап решения СЛАУ называется обратным ходом и состоит в последовательном определении неизвестных, начиная с x _n и заканчивая x ₁.

Точность результатов решения методом Гаусса определяется точностью арифметических операций при преобразовании элементов матрицы коэффициентов. Количество арифметических операций в методе Гаусса связано с размерностью системы и примерно равно (2/3)× n³. Контроль полученных решений производится путем их подстановки в исходную систему и вычисления невязок, т.е. разностей между правыми и левыми частями уравнений.

Благодаря ошибкам округления решение, получаемое с помощью точного прямого метода может сильно отличаться от истинного.

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

Для решения СЛАУ итерационными методами преобразуем систему уравнений от формы (8) к виду:

x ₁ = (b ₁ - a ₁₁ x ₁ - a ₁₂ x ₂ -... - a _1n x _n) / a ₁₁ + x ₁

x ₂ = (b ₂ - a ₂₁ x ₁ - a ₂₂ x ₂ -... - a _2n x _n) / a ₂₂ + x ₂ (9)

..................

x _n = (b _n - a _n1 x ₁ - a _n2 x ₂ -... - a _nn x _n) / a _nn + x _n

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

(10)

Процесс (10) можно видоизменить, если использовать приближения к решениям, найденные при выполнении текущей итерации. Как правило это приводит к ускорению сходимости. В этом случае итерационный метод называется методом Зейделя. Итерационный процесс имеет вид:

(11)

Для сходимости итерационных методов необходимо, чтобы значения диа-гональных элементов матрицы СЛАУ были преобладающими по абсолютной величине по сравнению с другими элементами. Условие сходимости можно обеспечить преобразованием матрицы коэффициентов путем перестановки уравнений и неизвестных.

Итерационный процесс заканчивают, когда выполнятся условия:

(12)

где e -заданная погрешность; k = 1, 2,..., n.

Ниже приведена программа решения системы уравнений методом Зейделя. В строках 20-50 программы вводятся значения числа уравнений N, погрешность вычислений E и максимально допустимое число итераций M. В строке 60 задаются начальные приближения переменных, равные нулю. В строках 70-100 вводятся коэффициенты A(I, J) и свободные члены B(I) системы уравнений. Итерационный цикл по переменной K реализуется в строках 170-230. Если за M итераций условия (12) не выполнятся, то выводится сообщение «Итерации все». Переменная L в случае, когда решение найдено, принимает значение счетчика итераций, в противном случае L=0. В строках 130-150 происходит вывод найденных значений переменных X(I) и числа итераций.

Для контроля и отладки программы решена система уравнений:

x ₁ + 4 x ₂+ 2 x ₃ = 3

2 x ₁ + 7 x ₂ + 3 x ₃ = 4

8 x₂ + x ₃ = -5

С абсолютной точностью 10^-5 за 4 итерации находится решение этой системы x ₁ = 1, x ₂ = -1, x ₃ = 3.

Программа решения системы линейных уравнений методом Зейделя:

10 PRINT " Решение системы линейных уравнений методом Зейделя."

20 INPUT " Задайте число уравнений N="; N

30 DIM A(N, N), B(N), X(N)

40 INPUT " Задайте погрешность вычислений E="; E

50 INPUT " Задайте максимально допустимое число итераций M="; M

60 FOR I=1 TO N: X(I)=0: NEXT I

70 FOR I=1 TO N

80 FOR J=1 TO N: PRINT " A" I; J;: INPUT A(I, J): NEXT J

90 PRINT " B" I;: INPUT B(I)

100 NEXT I

110 GOSUB 170

120 IF L=0 THEN 300

130 FOR I=1 TO N

140 PRINT " X" I" =" X(I): NEXT I

150 PRINT L" итераций"

160 GOTO 300

170 FOR K=1 TO M: L=K

180 FOR I=1 TO N: S=B(I)

190 FOR J=1 TO N: S=S-A(I, J)*X(J): NEXT J

200 S=S/A(I, I): X(I)=X(I)+S: IF ABS(S)> E THEN L=0

210 NEXT I

220 IF L< > 0 THEN RETURN

230 NEXT K

240 PRINT " итерации все"

250 RETURN

300 END

Расчет на ЭВМ:

Решение системы линейных уравнений методом Зейделя.

Задайте число уравнений N=? 3

Задайте погрешность вычислений E=? 0.00001

Задайте максимально допустимое число итераций M=? 20

A 1 1? 1

A 1 2? 4

A 1 3? 2

B 1? 3

X 1 = 1

X 2 = -1

X 3 = 3

4 итерации.

3. Интерполяция зависимостей