Решение нелинейных уравнений

Большинство уравнений описывающих реальные явления содержат нели­нейности. Под нелинейностями понимаются элементы, характеристики которых (y=F(x), где у - выход, x - вход) не являются линейными. На практике, как правило пользуются упрощенными моделями, например, пружинный маятник, нелинейность коэффициента упругости пружины считают постоянным, а это истинно только на небольшом интервале растяжений пружины. Но для нахождения точного решения этого не достаточно. Учет нелинейностей приводит к образованию сложных математических выражений, для которых не возможно записать решение в явной форме. В этом случае применяются приближенные, численные, методы решения.

Метод деления пополам

Метод основан на последовательном приближении к решению. Используются обозначения: FA=f(A), FB=f(B), FC=f(C).

1 -Выбирается отрезок на котором решение существует (FA>0)&&(FB<0) или (FA<0)&&(FB>0). Иначе решения нет.

2 - Отрезок [A B] делится пополам, получаем точку С.

3 - Проверяем является ли точка С решением с выбранной точностью (fabs(FC)<Точность). Если да, то решение найдено.

4 - Если знак FB совпадает с FA, как на рис., то принимаем за новое значение точки В точку С. Соответственно FB=FC. Возвращаемся к пп 2.

5 - Иначе перемещаем правую границу, те точка А. A=C, FA=FC. Переход к пп 2.

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

Метод секущих

Метод очень похож на предыдущий, иначе определяется только точка С. Можно предположить, что решение ближе к тому краю отрезка, где значение функции меньше. Проводим прямую линию [FA FB], и из двух подобных тре­угольников находим точку С. (|FA|/|FB|= (C-A)/(B-C). В остальном алго­ритм аналогичен предыдущему.

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

На практике часто встречаются уравнения, для которых нельзя ука­зать интервал (А-В) на котором находится решение. В током случае обыч­но используют метод Ньютона.

Из точки F проводится касательная к графику F(x), точка пересече­ния ее с осью x дает новое приближение решения X*. Это повторяется по­ка :

|F(x)| > <Точности решения>.

Угол наклона касательной можно найти через центральные разности:

При использовании данного метода необходимо учитывать возможность зацикливания и отсутствия решения. Эти случаи представленны ниже. Для того чтобы избежать эти ситуации рекомендуется ввести переменную intn_iterи увеличивать ее на каждой итерацииn_iter++. Если она превышает устновленный придел, то возможно решения нет.

Решение систем нелинейных уравнений

На практике задача отыскания решения системы уравнений встречает­ся чаще, чем решение уравнения с одним неизвестным. Рассмотрим для примера систему:

x1' + 10*(x1 - x2 )= 0

x2' + x1 * x3 -28*x1 +x2 = 0 или в векторной форме F(X)=0

x3' - x1 * x2 +2.6*x3= 0

Прямое решение такой системы невозможно. Единственный путь итерационные методы.

1.Выбирается некоторое приближение к решению (начальные значения)

2.Массиву предыдущих значений X присваиваем текущие значения X.

3.Вычисляется функция Якоби:

Эта функция является линеаризацией (упрощением) исходной системы уравнений в окрестностях точек x1,x2,x3

4. Считаем

F(X(n+1))=F(X(n)) + F'(X(n)) * (X(n+1)-X(n)) = 0

Или

F'(X(n)) * D X(n+1) = - F(X(n))

Решая эту систему относительно D X(n+1), получаем:

X(n+1) = X(n) + D X(n+1)

5. Полученные точки x1,x2,x3 подставляем в исходную систему уравне­ний. И если S | Fi(X) | <= Заданной точности, то решение найдено, ина­че возвращаемся к пп. 3.

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

7. Выводим или сохраняем рассчитанные точки, увеличиваем счетчик вре­мени. Если необходимо продолжать расчет, то возвращаемся к пп. 2.

Рекомендуется использовать:

float X[4] ; Значения независимых переменных *

float OldX[4] ; Значения переменных на предыдущем шаге *

float rX[4] ; Рабочий массив

float H ; Шаг интегрирования *

float dt ; Шаг дифференцирования *

float EPS ; Ошибка вычислений

float MaxEPS=0.05 ; Максимально допустимая ошибка вычислений

int N_ITER ; Счетчик количества итераций

int MaxITER=10 ; Максимально допустимое количество итераций

float A[4][4] ; Значения Якобиана

float B[4] ; Значения функции

float T ; Текущее время

float MaxT ; Общее время расчета

float G1[500],

G2[500],G3[500] ; Точки для графиков функции

float dX(int n)

{

return((X[n]-OldX[n])/H) ;

} /* Кон. dX() - Функция для вычисления производной n-ой перем. */

float F(int n)

{

float r ;

switch(n)

{

case 0:r=dX(0)+10*(X[0]-X[1]) ; break ;

case 1:r=dX(1)+X[0]*X[2]-28*X[0]+X[1] ; break ;

case 2:r=dX(2)+X[0]*X[1]+2.6*X[2] ; break ;

default:printf("\n Недопустимое значение функции %d",n) ;

} /* Кон. switch() */

} /* Кон. F() - Функция описания модели */

* - Рекомендуется описать как глобальные переменные.