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

max x_i⁽³⁾ - x_i⁽²⁾ .

¹ⁱ³

x₁⁽³⁾ - x₁⁽²⁾ = 2.152-1.583= 0.469

x₂⁽³⁾ - x₂⁽²⁾ = - 1.945+ 1.569= 0.376

x₃⁽³⁾ - x₃⁽²⁾ = 2.024-1.788= 0.236

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

Требуется решить систему нелинейных уравнений вида:

F₁(x₁,x₂,x₃,...,x_n) = 0

F₂(x₁,x₂,x₃,...,x_n) = 0

. . . . . . . . . . . (3.1)

F_n(x₁,x₂,x₃,...,x_n) = 0

3.1 Метод простой итерации (метод Якоби)

Систему нелинейных уравнений (3.1) после преобразований

x_i =x_i - F_i(x)/M_i ; i=1, 2, 3, ..., n

(здесь M_i определяются из условия сходимости), представим в виде:

x₁ =f₁(x₁,x₂,x₃,...,x_n)

x₂ =f₂(x₁,x₂,x₃,...,x_n)

. . . . . . . . . . . (3.2)

x_n =f_n(x₁,x₂,x₃,...,x_n)

Из системы (2.3) легко получить итерационные формулы метода Якоби. Возьмем в качестве начального приближения какую-нибудь совокупность чисел x₁⁽⁰⁾,x₂⁽⁰⁾,...,x_n⁽⁰⁾. Подставляя их в правую часть (3.2) вместо переменных x₁,x₂,...,x_n получим новое приближение к решению исходной системы:

x₁⁽¹⁾ =f₁(x₁⁽⁰⁾,x₂⁽⁰⁾,x₃⁽⁰⁾,...,x_n⁽⁰⁾)

x₂⁽¹⁾ =f₂(x₁⁽⁰⁾,x₂⁽⁰⁾,x₃⁽⁰⁾,...,x_n⁽⁰⁾)

. . . . . . . . . . . (3.3)

x_n⁽¹⁾ =f_n(x₁⁽⁰⁾,x₂⁽⁰⁾,x₃⁽⁰⁾,...,x_n⁽⁰⁾)

Эта операция получения первого приближения x₁⁽¹⁾,x₂⁽¹⁾,...,x_n⁽¹⁾ решения системы уравнения (3.2) называется первым шагом итерации. Подставляя полученное решение в правую часть уравнения (3.2) получим следующее итерационное приближение: (x₁⁽²⁾,x₂⁽²⁾,...,x_n⁽²⁾) и т.д.

Итерационный процесс можно считать законченным, если все значения переменных, полученных k+1-ой итерации, отличается от значений соответствующих переменных, полученных от предыдущей итерации, по модулю меньше наперед заданной точности , т.е. если:

maxx_i^(k+1) -x_i^(k) . (3.4)

¹ⁱⁿ

3.2 Метод Зейделя

Метод Зейделя отличается от метода Якоби тем, что вычисления ведутся не по формулам (3.3), а по следующим формулам:

x₁^(k+1) =f₁(x₁^(k),x₂^(k),x₃^(k),....,x_n^(k))

x₂^(k+1) =f₂(x₁^(k+1),x₂^(k),x₃^(k),....,x_n^(k))

. . . . . . . . . . . (3.5)

x_n^(k+1) =f_n(x₁^(k+1),x₂^(k+1),x₃^(k+1),...,x_n^(k)).

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

Сходимость метода Зейделя (Якоби тоже) зависит от вида функции в (3.2), вернее она зависит от матрицы составленной из частных производных:

f₁₁ f₁₂ f₁₃ ... f_1n

F = f₂₁ f₂₂ f₂₃ ... f_2n

. . . . . . . . . . . (3.6)

f_n1 f_n2 f_n3 ... f_nn

где f_ij = f_i(x₁,x₂,...x_n)/x_j.

Итерационный процесс сходится, если сумма модулей каждой строки F меньше единицы в некоторой окрестности корня:

f_i1f_i2f_i3....f_in ; i=1, 2, 3, ..., n

_n

или max f_ij .

^{1in j =1}

Пример: Найти решение системы (3.7) методом Зейделя с точностью =0,001.

F(x,y)=2sin(x+1)-y-0.5 = 0 (3.7)

G(x,y)=10cos(y-1)-x+0.4 = 0

Решение: Представим (3.7) в виде (3.5):

x = f₁(x,y) = x-(2sin(x+1)-y-0,5)/M₁ (3.8)

y = f₂(x,y) = y-(10cos(y-1)-x+0,4)/M₂

Строим графики кривых системы (3.7) и определяем начальные приближения x₀=-1, y₀=-0,7.

Запишем достаточное условие сходимости и определяем M₁, M₂:

f_1x f_1y 1-2cos(x+1)/ M₁ 1/ M₁

F = =

f_2x f_2y -1/M₂ 1+10sin(y-1)/M₂

1-2cos(x₀+1)/M₁+1/M₁ и -1/M₂+1+10sin(y₀-1)/M₂ ,

1-2cos(1-1)/M₁+1/M₁ и -1/M₂+1+10sin(-0,7-1)/M₂ ,

1-2/M₁+1/M₁ и -1/M₂+ 1-9,91665/M₂ ,

Определяем частные значения M₁=2, M₂=10, которые удовлетворяют неравенствам

1-2/2+1/2 и 1/10+ 1-9,91665/10 ,

Переходим к реализации итерационного процесса:

x_k+1 = x_k-(2·sin(x_k+1)-y_k-0,5)/2

y_k+1 = y_k-(10·cos(y_k -1)-x_k+1+0,4)/10

x₁=x₀-(2·sin(x₀+1)-y₀-0,5)/2= -1-(2·sin(-1+1)+0,7-0,5)/2=-1,1

y₁=y₀-(10·cos(y₀-1)-x₁+0,4)/10=

=-0,7-(10·cos(-0,7-1)+ 1,1+0,4)/10=-0,72116

x₂=x₁-(2·sin(x₁+1)-y₁-0,5)/2=

=-1,1-(2·sin(-1,1+1)+ 0,72116-0,5)/2=-1,11075

y₂=y₁-(10·cos(y₁-1)-x₂+0,4)/10=

=-0,72116-(10·cos(-0,72116-1)+1,11075+0,4)/10=-0,72244

x₃=x₂-(2·sin(x₂+1)-y₂-0,5)/2=

=-1,11075-(2·sin(-1,11075+1)+ 0,72244-0,5)/2=-1,11145

y₃=y₂-(10·cos(y₂-1)-x₃+0,4)/10=

=-0,72244-(10*cos(-0,72244-1)+1,11145+0,4)/10=-0,72252

Определяем погрешность по формуле

maxx_i^(k+1) -x_i^(k) .

¹ⁱⁿ

x₃ -x₂ -1,11145+1,11075= 0,0007 =0,001

y₃ - y ₂ -0,72252+0,72244= 8E-05 =0,001

Таким образом, имеем решение: x =-1,11145; y =-0,72252.

Программа, реализующая решение данной задачи, представлена на рис.3.2.

CLS

REM LR-3-2, m=13, n=5

INPUT X,Y, M1,M2

1 X=X-(2*SIN(X+1)-Y - 0.5)/M1

Y=Y-(10*COS(Y-1)-X+0.4)/M2

PRINT X,Y

INPUT TT

GOTO 1

END

Рис.3.2. Программа решения методом Зейделя.

3.3 Метод Ньютона

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

F(x, y)=0

G(x, y)=0 (3.9)

Пусть известно некоторое приближение x_k,y_k корня x,y. Тогда поправки x_k=x_k+1-x_k, y_k+y_k+1-y_k можно найти решая систему:

F(x_k+x_k, y_k+y_k)=0

G(x_k+x_k, y_k+y_k)=0. (3.10)

Для этого разложим функции F, G в ряд Тейлора по x_k, y_k. Сохранив только линейные по x_k, y_k части, получим систему линейных уравнений:

F(x_k,y_k) F(x_k,y_k)

x_k + y_k = - F(x_k,y_k)

x y

(3.11)

G(x_k,y_k) G(x_k,y_k)

x_k + y_k = - G(x_k,y_k)

x y

относительно неизвестных поправок x_k и y_k. Решая эту систему линейных уравнений, определяем значения x_k, y_k.

Таким образом, решение системы уравнений по методу Ньютона состоит в построении итерационной последовательности:

x_k+1 =x_k +x_k

y_k+1 =y_k +y_k (3.12)

где x_k, y_k - решения систем линейных уравнений, вида (3.11) на каждом шаге итерации. В методе Ньютона также для обеспечения хорошей сходимости важен правильный выбор начального приближения.

Пример: Найти решение системы (3.7) методом Ньютона с точностью =0,001.

F(x,y)=2·sin(x+1)-y-0.5 = 0

G(x,y)=10·cos(y-1)-x+0.4 = 0 (3.13)

методом Ньютона, при этом добиться точности 0.001.

Решение: Начальные приближения x₀=-1 и y₀=-0,7. Определим частные производные:

F(x, y) F(x, y)

= 2·cos(x+1); = -1;

x y

G(x, y) G(x, y)

-1; = -10·sin(y-1).

x y

и используя (3.11), построим систему линейных уравнений относительно поправок

2·cos(x₀+1) · x₀ -1 · y₀ = - 2·sin(x₀+1)+y₀+0.5

-1· x₀ -10·sin(y₀ -1) · y₀ = - 10·cos(y₀ -1)+x₀ -0.4

Подставляя начальные приближения x₀=-1; y₀=-0,7 и решая систему линейных уравнений определяем поправки на первом шаге итерации

x₀=-0,1112 y₀=-0,0225

Далее начальное приближение уточняем по формулам (3.12)

x₁ =x₀ +x₀= - 1 - 0,1112 = - 1,1112

y₁ =y₀ +y₀= - 0,7 - 0,0225= - 0,7225

Подставляя результаты первой итерации x₁=-1,1112; y₁=-0,7225 и решая систему линейных уравнений

2·cos(x₁+1) · x₁ -1 · y₁ = - 2·sin(x₁+1)+y₁+0.5

-1· x₁ -10·sin(y₁ -1) · y₁ = - 10·cos(y₁ -1)+x₁ -0.4

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

x₁=0,00026 y₁=0,00004

Далее x₁ и y₁ уточняем по формулам (3.12)

x₂ =x₁ +x₁= -1,1112+0,00026= - 1,1115

y₂ =y₁ +y₁=-0,7225+0,00004= - 0,7225

Определяем погрешность по формуле

maxx_i^(k+1) -x_i^(k) .

¹ⁱⁿ

x₂ -x₁ или x₁ 0,00026 =0,001

y₂ - y₁ или y₁ 0,00004 =0,001

Таким образом имеем решение: x = - 1,1115; y = - 0,7225.

Программа реализующая метод Ньютона для указанной задачи представлена на рис. 3.3.

REM LR-3-3, m=13, n=5

INPUT X, Y

1 F = 2*SIN(X+1)-Y - 0.5

G = 10*COS(Y-1)-X+0.4

Fx =2*COS(X+1)

Fy =-1

Gx =-1

Gy =-10*SIN(Y-1)

D = Fx*Gy- Gx*Fy

DX=(G*Fy-F*Gy)/D

DY=(F*Gx-G*Fx)/D

X =X+DX

Y =Y+DY

PRINT X;Y; F;G;DX;DY;

INPUT TT

GOTO 1

END

Рис.3.3. Программа, реализующая метод Ньютона

4.1 Приближение функции по методу наименьших квадратов (МНК)

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

Как правило, общий вид этой функциональной зависимости или так называемой эмпирической формулы известен, а некоторые числовые параметры закона неизвестны.

Процесс выражения опытных данных функциональной зависимостью с помощью метода наименьших квадратов состоит из двух этапов: на первом этапе выбирают вид искомой формулы, а на втором этапе для формулы подбирают параметры. Для первого этапа удобно графическое представление этой зависимости и по ней выбрать вид возможной зависимости во втором этапе, в соответствии с идеей МНК, необходимо минимизировать сумму отклонений:

S = (y(x_i)-y_i)² (4.1)

где x_i, y_i - значения опытных данных y(x_i) - значение функции, вычисленное в точке x_i ; n - число опытов.

В случае линейной эмпирической формулы (4.1) принимает вид:

S(a,b) = (ax_i + b - y_i)² min (4.2)

Функция (4.2) имеет минимум в точках, в которых частные производные от S по параметрам a и b обращаются в нуль, т.е.

S(a,b)/a=0, S(a,b)/b=0 (4.3)

2(a x_i+ b - y_i) x_i = 0

2(a x_i+ b - y_i) = 0

ax_i² + bx_i =x_i y_i

ax_i + b n = y_i (4.4)

Решая систему уравнений (4.4), получим значения a и b уравнения y=ax+b. Пример: Подобрать аппроксимирующий полином первой степени y=ax+b для данных

x_i 0 1 2 3

y_i 0.1 0.9 2.1 3

Решение: Для удобства вычисленные значения расположим в таблице.

Система для определения коэффициентов имеет вид:

14 a + 6 b = 14.1

6 a + 4 b = 5.1 (4.5)

Решая систему (4.5), получим следующие значения параметров a=1.29; b=-0.675. Следовательно, искомый полином имеет вид y= 1.29 x - 0.67.

В случае квадратичной зависимости (4.1) принимает вид:

S(a,b,c) = (ax_i² + bx_i + c - y_i)² min (4.6)

Функция (4.6) имеет минимум в точках, в которых частные производные от S по параметрам a, b, c обращаются в нуль:

S(a,b,c)/a=0, S(a,b,c)/b=0, S(a,b,c)/c=0 (4.7)

В результате дифференцирования и элементарных преобразований для определения параметров получают систему линейных уравнений:

S(a,b,c)/a = 2 (ax_i² + bx_i + c - y_i) x_i²=0

S(a,b,c)/b = 2 (ax_i² + bx_i + c - y_i) x_i=0

S(a,b,c)/c = 2 (ax_i² + bx_i + c - y_i) *1=0

Система линейных уравнений состоит из трех уравнений с тремя неизвестными:

ax_i⁴ + bx_i ³ + c x_i² =x_i² y_i

ax_i³ + bx_i ² + c x_i =x_i y_i (4.8)

ax_i² + bx_i + c n = y_i

Решая систему линейных уравнений (4.8) получим значения a, b, c уравнения y=ax²+bx+c.

Пример: Используя МНК построить эмпирическую зависимость y=ax²+bx+c, аппроксимирующую следующие табличные значения:

x_i -2 -1 0 1 2

y_i 6 2 -1 -2 -1

Тогда система линейных уравнений (4.8) относительно значений a, b, c примет вид:

34 a+ 0 b+10 c= 20

0 a+10 b+ 0 c= -18 (4.9)

10 a+ 0 b+ 5 c = 4

Решая систему (4.9) получим следующие значения параметров a=0.857; b=-1.8; c=-0.914. Таким образом, искомый полином имеет вид:

y=0.857 x² - 1.8 x - 0.914.

По таблице 4.3 можно определить сумму квадратов отклонений экспериментальных данных y_i от расчетных y_i S=0.112.

Каким образом можно уменьшить погрешность S ? Для этого необходимо увеличить число коэффициентов полинома y_. Это приведет к увеличению размерности системы. Таким образом, стараясь улучшить расчеты, т.е. уменьшить погрешность, приходится увеличивать объем вычислений.

4.2 Интерполяционный полином в форме Лагранжа

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

Пусть на некотором промежутке [a,b] заданы n различных узлов x₁,x₂,x₃, ..., x_n, а также значения некоторой функции y₁,y₂,y₃, ..., y_n в этих узлах. Необходимо построить полином P(x), проходящий через заданные точки, т.е.

P(x_i)=y_i

Этот полином называется интерполяционным полиномом, является единственным полиномом степени не выше n-1, и может быть записан, например, в форме Лагранжа или Ньютона.

Интерполяционный полином Лагранжа имеет следующую формулу:

P(x)= L_n-1(x)=_i l_i(x), (4.10)

(x-x₁)...(x-x_i-1)(x-x_i+1)...(x-x_n)

где l_i(x) =

(x_i-x₁)...(x_i-x_i-1)(x_i-x_i+1)...(x_i-x_n)

фундаментальные полиномы Лагранжа. Они удовлетворяют равенствам

1, если i = k

l_k(x_i) = (4.11)

0, если i k ,

и зависят лишь от заданных узлов x_i , но не от интерполируемой функции y_i .

Необходимо построить интерполяционный полином Лагранжа, проходящий через заданные точки

L_n-1(x_i)=y_i , i=1,2,...,n.

Решение. Запишем фундаментальные полиномы Лагранжа:

(x-x₂)(x-x₃)(x-x₄) (x-0)(x-1/2)(x-1)

l₁(x)= = = - (x³-3/2 x²+1/2 x)/3,

(x₁-x₂)(x₁-x₃)(x₁-x₄) (-1-0)(-1-1/2)(-1-1)

(x-x₁)(x-x₃)(x-x₄) (x+1)(x-1/2)(x-1)

l₂(x)= = = 2(x³-1/2 x²- x+1/2),

(x₂-x₁)(x₂-x₃)(x₂-x₄) (0+1)(0-1/2)(0-1)

(x-x₁)(x-x₂)(x-x₄) (x+1)(x-0)(x-1)

l₃(x)= = = -8(x³- x)/3,

(x₃-x₁)(x₃-x₂)(x₃-x₄) (1/2+1)(1/2-0)(1/2-1)

(x-x₁)(x-x₂)(x-x₃) (x+1)(x-0)(x-1/2)

l₄(x)= = = 1(x³+1/2 x²- 1/2x).

(x₄-x₁)(x₄-x₂)(x₄-x₃) (1+1)(1-0)(1-1/2)

Например, для l₄(x) можно проверить свойство (5.2).

l₄(-1)=0, l₄(0)=0, l₄(1.2)=0, l₄(1)=1.

Подставляя l_i(x) в полином Лагранжа находим:

L₃(x)=y₁ l₁(x)+y₂ l₂(x)+y₃ l₃(x)+y₄ l₄(x) =

= 0 l₁(x)+2 l₂(x)+9/8 l₃(x)+0 l₄(x) = x³ - 2 x² - x + 2.

4.3 Интерполяционный полином в форме Ньютона

Запишем интерполяционный полином в форме Ньютона:

N_n-1(x)= y₁+¹(x₁,x₂)(x-x₁)+ ²(x₁,x₂,x₃)(x-x₁)(x-x₂)+...

+ ^n-1(x₁,x₂,...,x_n-1)(x-x₁)...(x-x_n-1), где

⁰_i = y_i,

¹(x_i,x_k)=(⁰_i - ⁰_k)/(x_i-x_k) - разделенная разность первого порядка,

²(x_i,x_j,x_k)=(¹(x_i,x_j)-¹(x_j,x_k))/(x_i-x_k) - разделенная разность второго порядка,

³(x_i,x_j,x_l,x_k)=(²(x_i,x_j,x_l)-²(x_j,x_l,x_k))/(x_i-x_k) - разделенная разность третьего порядка и т.д.

Для заданной таблицы 4.4 построим интерполяционный полином в форме Ньютона

N₃(x_i)=y_i .

Расчеты представим в виде таблицы 4.5.

Таблицы 4.5

i	xi	yi	1(i,k)	2(i,j,k)	3(i,j,l,k)
1	-1	0
			2
2	0	2		-5/2
			-7/4		1
3	1/2	9/8		-1/2
			-9/4
4	1	0

¹(1,2)=(y₁ - y₂)/(x₁ -x₂)=(0-2)/(-1-0)=2,

¹(2,3)=(y₂ - y₃)/(x₂ -x₃)=(2-9/8)/(0-1/2)= -7/4,

¹(3,4)=(y₃ - y₄)/(x₃ -x₄)=(9/8-0)/(1/2-1)= -9/4,

²(1,2,3)=(¹(1,2)-¹(2,3))/(x₁ -x₃)=(2+7/4)/(-1-1/2)= -5/2,

²(2,3,4)=(¹(2,3)-¹(3,4))/(x₂ -x₄)=(-7/4+9/4)/(0-1)= -1/2,

³(1,2,3,4)=(²(1,2,3)-²(2,3,4))/(x₁-x₄)=(-5/2+1/2)/(-1-1)=1,

N₃(x)=(1)+(1,2)(x-x₁)+(1,2,3)(x-x₁)(x-x₂)+

+(1,2,3,4)(x-x₁)(x-x₂)(x-x₃)=

=0+2(x+1)-5/2(x+1)x+1(x+1)(x-1/2)=x³ -2x²-x+2.

ЛИТЕРАТУРА

1. Калиткин Н.П. Численные методы. М.: Наука, 1978. - 512 с.

2. Солодовников А.С. Введение в линейную алгебру и линейное программирование. М.: Просвещение, 1966. - 183 с.

3. Химмельблау Д. Прикладное нелинейное программирование. М.: Мир, 1975. - 534 с.

4. Попов В.И. Численные методы расчета мостовых конструкций на ЭВМ. М.: 1981. - 78 с.

5. Монахов В.М. и др. Методы оптимизации. Применение математических методов в экономике. М., Просвещение, 1978. - 175 с.

6. Информатика. Методические указания к лабораторным работам. «Численные методы решения задач строительства на ЭВМ» для студентов специальности 2910.//Казанская государственная архитектурно-строительная академия; Сост.: Габбасов Ф.Г., Гатауллин И.Н., Киямов Х.Г. - Казань:, 1998. -50 с.

Размещено на Allbest.ru