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

ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ ЗАДАЧ

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

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

Этап 1. Исключение переменной . Разделим первое уравнение системы на первый коэффициент 2.3, получим:

отсюда выражаем .

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

отсюда систем имеет вид:

тогда переменная исключена.

Этап 2. Теперь исключим переменную . Разделим уравнение 2 на коэффициент перед , получим:

после подсчета коэффициентов получим:

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

отсюда получаем систему:

Далее выполняем последовательно, начиная с третьего и заканчивая первым, вычисление каждой переменой с помощью данных уравнений:

После вычисления переменных, мы получили решение системы с помощью метода Гаусса.

Общая методика метода Гаусса.

Система представлена в общем виде как:

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

Рассмотрим этапы по шагам.

Среди элементов , i, j=1,..., n матрицы A выбирается наибольший по модулю , называемый главным элементом. Соответствующая строка матрицы A с номером p называется главной строкой. Затем меняется местами первая строка со строкой p и первый столбец со столбцом q и осуществляется перенумерация коэффициентов и неизвестных. Информация о перенумерации запоминается. В результате этих операций первая строка становится главной.

Далее первое уравнение системы делят на = :

где .

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

где , , i=2,..., n; j=2,..., n

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

Объединяют все главные уравнения, в результате чего образуется система уравнений с треугольной матрицей, эквивалентной исходной:

После этого проводится обратный ход: из системы неизвестные определяются в обратном порядке по формулам:

…

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

Пример. Написать алгоритм и программу для решения системы линейных уравнений

методом Гаусса.

Программа выполняется следующим образом:

1. Из файла com.dat считываем коэффициенты системы в массив для решения системы. Данные в файле даны в следующем виде:

2.3 4.5 3.1 8

0.2 2.2 5 1

4 -2 7.8 -5

2. Проводим вычисления с помощью процедуры gauss, выполняющую решение системы линейных уравнений методом Гаусса с выбором главного элемента. Для ее адекватной работы в процедуру следует передать переменные mas типа tablica (размеры задаются по потребности в разделе типов), в которой находятся коэффициенты системы линейных уравненений. Вторая целочисленная переменая cislo передает размерность системы (в нашем случае 3). После служебного слова var передаются переменные значения которых возвращает процедура, т.е. она изменяет их значения. В переменной X типа colonka (размеры должны соответствовать размерности системы и задаются в разделе определения типов) после выполнения процедуры будет находится решение системы линейных уравнений (по порядку будут находится значения переменных). Переменная flag_virogdena типа boolean показывает если она true, то система линейных уравнений является вырожденной, ели она false, то система линейных уравнений не вырождена.

3. После решения системы линейных уравнений с помощью процедуры gauss, выводим решения системы на экран.

uses crt;

type

tablica = array[0..2, 0..3] of real;

colonka = array[0..2] of real;

var

cislo: integer;

temp: array[0..3] of real;

znaki: array[0..2] of integer;

f: text;

A: tablica;

X: colonka;

flag: boolean;

i, j: integer;

procedure gauss(mas: tablica; cislo: integer; var X: colonka; var flag_virogdena: boolean);

var i, j, k, index, ind, znak: integer;

sum, max: real;

begin {начало процедуры}

flag_virogdena: =false;

for i: =0 to cislo-1 do

znaki[i]: =i;

for k: =1 to cislo-1 do begin

max: =mas[k-1, k-1];

index: =k-1; ind: =k-1;

for i: =k-1 to cislo-1 do

for j: =k-1 to cislo-1 do

if abs(mas[i, j])> abs(max) then begin

max: =mas[i, j];

index: =i;

ind: =j

end;

if max=0 then flag_virogdena: =true;

if not(flag_virogdena) then begin

for j: =k-1 to cislo-1 do begin

temp[j]: =mas[index, j];

mas[index, j]: =mas[k-1, j];

mas[k-1, j]: =temp[j];

end;

temp[cislo]: =mas[index, cislo];

mas[index, cislo]: =mas[k-1, cislo];

mas[k-1, cislo]: =temp[cislo];

for i: =0 to cislo-1 do begin

temp[i]: =mas[i, ind];

mas[i, ind]: =mas[i, k-1];

mas[i, k-1]: =temp[i];

end;

znak: =znaki[ind];

znaki[ind]: =znaki[k-1];

znaki[k-1]: =znak;

for i: =k to cislo-1 do

mas[k-1, i]: =mas[k-1, i]/mas[k-1, k-1];

mas[k-1, cislo]: =mas[k-1, cislo]/mas[k-1, k-1];

mas[k-1, k-1]: =1;

for i: =k to cislo-1 do begin

for j: =k to cislo-1 do

mas[i, j]: =mas[i, j]-mas[i, k-1]*mas[k-1, j];

mas[i, cislo]: =mas[i, cislo]-mas[i, k-1]*mas[k-1, cislo];

mas[i, k-1]: =0;

end;

end;

end;

for i: =0 to cislo-1 do

X[i]: =0;

for k: =0 to cislo-1 do begin

sum: =0;

if (k=0) and (mas[cislo-1-k, cislo-1-k]=0) then flag_virogdena: =true;

if (cislo-1-k)< cislo-1 then

for j: =cislo-k to cislo-1 do

sum: =sum-mas[cislo-1-k, j]*X[j];

if mas[cislo-1-k, cislo-1-k]=0 then X[cislo-1-k]: =0 else

X[cislo-1-k]: =(mas[cislo-1-k, cislo]+sum)/mas[cislo-1-k, cislo-1-k];

end;

for i: =0 to cislo-1 do

temp[i]: =X[i];

for i: =0 to cislo-1 do

X[znaki[i]]: =temp[i];

end; {конец программы}

Begin{н.п} {программа для решения системы уравнений}

assign(f, 'com.dat');

reset(f);

for i: =0 to 2 do begin{н.б.1}

for j: =0 to 3 do

read(f, a[i, j]); {Ввод коэффициентов системы}

readln(f);

end; {к.б.1}

gauss(a, 3, x, flag); {вычисление решения }

for i: =0 to 2 do

writeln('x', i+1, '=', x[i]: 7: 2);

close(f);

readln;

End.{к.п.}

Результаты задачи 1:

Файл com.dat:

2.3 4.5 3.1 8

0.2 2.2 5 1

4 -2 7.8 -5

x1=0.83

x2=1.77

x3=-0.61

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

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

(4.2.1)

где x- неизвестная величина; - параметры

Для численного решения уравнения (4.2.1) одним из эффективных является метод дихотомии. Сущность этого метода состоит в поэтапном делении отрезка [a, b] попалам, на котором находится решение уравнения (4.2.1).

На первом этапе табулируют функцию на отрезке с точностью до выбранного шага решения уравнения (4.2.1) при заданных параметрах . В результате получают таблицу:

				...
				...

Из таблицы выделяют интервалы на которых функция изменяет свой знак (F(x)> 0 и F(x)< 0), что проиллюстрировано на рисунке

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

Метод дихотомии. Рассмотрим идеалогию метода дихотомии. Первоначально делят отрезок [a, b] на котором функция F(x) изменяет свой знак попалам, т.е. выбирают начальное приближение корня равным . Если F()=0, то является корнем уравнения. В противном случае анализируют тот из отрезков [a, ] или [ , b], на концах которого функция f(x) имеет разные знаки. Далее полученный отрезок вновь делят пополам и действия повторяют до тех пор пока не будет найден точно корень уравнения либо не будет достигнута заданная точность (т.е. пока длина отрезка не станет меньше 2e). В этом случае середина последнего отрезка дает значение корня с требуемой точностью.

Иллюстрация метода.

1-й шаг. Вычисление и проверка условия F(a)*F()=0: если оно не верно, тогда проверяется условие F(a)*F()< 0. Если оно верно тогда рассматривается отрезок [a, ] и проверяется критерий остановки поиска корня . При не удовлетворении условия выполняется шаг 2.

2-й шаг. Вычисление (деление отрезка попалам). Проверка условия F(a)*F()=0. Если оно не верно, тогда проверяется условие F(a)*F()< 0. Если оно верно, тогда рассматривается отрезок [a, ] и проверяется критерий остановки . При не удовлетворении условия выполняется шаг 3.

3-й шаг. Вычисление . Проверка условия F(a)*F()=0. Если оно верно, значит -корень уравнения.

Подробная блок–схема алгоритма дихотомии представлена в примере данного подраздела.

Метод Ньютона. Для использования данного метода F(x)-должна быть дифференцируема. Задается начальное значение переменной , причем необходимо выбирать из условия . Поиск значения x ведется по формуле . При выполнении условия окончания вычисления заканчиваются, где q-коэффициент наперед заданный.

Пример 1. Составить схему алгоритма метода Ньютона.

Обозначим в схеме как “условие 1”.

Пример 2. Написать алгоритм и программу решения нелинейного уравнения на отрезке [-10, 4] методом деления попалам. Построить график. Выбрать отрезок для уточнения корня.

Для вычисления пределов изменения корней протабулируем уравнение F(x)=0 на отрезке [-10, 4](выбор интервала произволен).

Из таблицы видно, что на отрезках [-10, -9] [-7, -6] [-4, -3] [0, 1] [2, 3] меняется знак функции, значит на этих отрезках есть корни уравнения. При состовлении программы метода дихотомии найдем решение уравнения на отрезке [0, 1].

График функции будет иметь вид:

Обозначим в схеме алгоритма выражения

как “Условие1”, как “Условие 2”,

как “Условие3”

program noline_eqution; {решение нелинейного уравнения}

uses crt;

var

xn, xk, epsilon, x: real;

flag_stop: boolean;

function fun(x: real): real;

begin

fun: =exp(x)-10*sin(x);

end;

begin {н.п.}

clrscr;

xn: =0;

xk: =1;

epsilon: =0.005;

flag_stop: =false;

if fun(xn)*fun(xk)> 0 then writeln('Решения нет на данном отрезке')

else begin

while (abs(xk-xn)> 2*epsilon) and (not(flag_stop)) do begin

x: =(xn+xk)/2;

if fun(x)=0 then flag_stop: =true

else if fun(xn)*fun(x)< 0 then xk: =x else xn: =x;

writeln('x=', x: 5: 4, ’ xn=’, x0: 3: 2, ’ xk=’, xn: 3: 2);

end;

end;

readln;

end. {к.п.}

Результаты задачи 2:

x=0.5000 xn=0.00 xk=0.50

x=0.2500 xn=0.00 xk=0.25

x=0.1250 xn=0.00 xk=0.13

x=0.0625 xn=0.06 xk=0.13

x=0.0938 xn=0.09 xk=0.13

x=0.1094 xn=0.11 xk=0.13

x=0.1172 xn=0.11 xk=0.12