Численные методы решения типовых математических задач

1.1 Постановка задачи

Разработать схему алгоритма и написать программу на языке Turbo Pascal 7.0 для решении систем линейных алгебраических уравнений, используя метод простой итерации.

1.2 Математическая формулировка задачи

Пусть А - невырожденная матрица и нужно решить систему

где диагональные элементы матрицы А ненулевые.

1.3 Обзор существующих численных методов решения задачи

Метод Гаусса

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

Пусть дана СЛАУ

Запишем расширенную матрицу системы:

На первом шаге алгоритма Гаусса выберем диагональный элемент (если он равен 0, то первую строку переставляем с какой-либо нижележащей строкой) и объявляем его a₁₁?0 ведущим, а соответствующую строку и столбец, на пересечении которых он стоит - ведущими. Обнулим элементы ведущего столбца. Для этого сформируем числа (-a₂₂/a₁₁), (-a₃₁/a₁₁), .. , (a_n1/a₁₁).

LU-разложения матриц

Компьютерная реализация метода Гаусса часто осуществляется с использованием LU-разложения матриц.

LU - разложение матрицы A представляет собой разложение матрицы A в произведение нижней и верхней треугольных матриц, т.е.

A=LU

где L - нижняя треугольная матрица (матрица, у которой все элементы, находящиеся выше главной диагонали равны нулю, lij=0 при i>j), U- верхняя треугольная матрица (матрица, у которой все элементы, находящиеся ниже главной диагонали равны нулю, uij=0 при i>j).

LU - разложение может быть построено с использованием описанного выше метода Гаусса. Рассмотрим k - ый шаг метода Гаусса, на котором осуществляется обнуление поддиагональных элементов k - го столбца матрицы. Как было описано выше, с этой целью используется следующая операция

В терминах матричных операций такая операция эквивалентна умножению A(k)=MkA(k-1), где элементы матрицы определяются следующим образом

В результате прямого хода метода Гаусса получим , A(n-1)=U

где A(n-1)=U - верхняя треугольная матрица, а - нижняя треугольная матрица, имеющая вид .

Таким образом, искомое разложение A=LU получено.

Метод прогонки

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

решение которой будем искать в виде

где Qi,Pi,i=1,n - прогоночные коэффициенты, подлежащие определению. Для их определения выразим из первого уравнения СЛАУ (1.1) x1 через x2, получим:

откуда

Из второго уравнения СЛАУ (1.1) с помощью (1.3) выразим x2 через x3, получим:

откуда

Продолжая этот процесс, получим из i-го уравнения СЛАУ (1.1):

следовательно

Из последнего уравнения СЛАУ имеем

то есть

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

При большом числе уравнений прямые методы решения СЛАУ (за исключением метода прогонки) становятся труднореализуемыми на ЭВМ прежде всего из-за сложности хранения и обработки матриц большой размерности. В то же время характерной особенностью ряда часто встречающихся в прикладных задачах СЛАУ является разреженность матриц. Число ненулевых элементов таких матриц мало по сравнению с их размерностью. Для решения СЛАУ с разреженными матрицами предпочтительнее использовать итерационные методы.

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

Метод Зейделя решения СЛАУ

Метод простых итераций довольно медленно сходится. Для его ускорения существует метод Зейделя, заключающийся в том, что при вычислении компонента xik+1 вектора неизвестных на (k+1)-й итерации используются x1k+1, x2k+1, .., xik+1 уже вычисленные на (k+1)-й итерации. Значения остальных компонент берутся из предыдущей итерации. Так же как и в методе простых итераций строится эквивалентная СЛАУ и за начальное приближение принимается вектор правых частей . Тогда метод Зейделя для известного вектора на k-ой итерации имеет вид:

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

Из этой системы видно, что , где В - нижняя треугольная матрица с диагональными элементами , равными нулю, а С - верхняя треугольная матрица с диагональными элементами, отличными от нуля, б=B+C. Следовательно

При таком способе приведения исходной СЛАУ к эквивалентному виду метод простых итераций носит название метода Якоби.

откуда

Таким образом, метод Зейделя является методом простых итераций с матрицей правых частей б=(E-B)-1C и вектором правых частей (E-B)-1в.

Разрешим систему относительно неизвестных при ненулевых диагональных элементах , aii?0, i= 1,n (если какой-либо коэффициент на главной диагонали равен нулю, достаточно соответствующее уравнение поменять местами с любым другим уравнением). Получим следующие выражения для компонентов вектора в и матрицы б эквивалентной системы

В качестве нулевого приближения x⁽⁰⁾вектора неизвестных примем вектор правых частей x⁽⁰⁾=в. Тогда метод простых итераций примет вид:

Из (1.19) видно преимущество итерационных методов по сравнению, например, с рассмотренным выше методом Гаусса. В вычислительном процессе участвуют только произведения матрицы на вектор, что позволяет работать только с ненулевыми элементами матрицы, значительно упрощая процесс хранения и обработки матриц.

Имеет место следующее достаточное условие сходимости метода простых итераций.

Метод простых итераций (1.19) сходится к единственному решению СЛАУ при любом начальном приближении x⁽⁰⁾, если какая-либо норма матрицы б эквивалентной системы меньше единицы

Если используется метод Якоби (выражения (1.18) для эквивалентной СЛАУ), то

достаточным условием сходимости является диагональное преобладание матрицы A, т.е.

(для каждой строки матрицы A модули элементов, стоящих на главной диагонали, больше суммы модулей недиагональных элементов). Очевидно, что в этом случае ||б||_c меньше единицы и, следовательно, итерационный процесс (1.19) сходится.

Приведем также необходимое и достаточное условие сходимости метода простых итераций. Для сходимости итерационного процесса (1.19) необходимо и достаточно, чтобы спектр матрицы б эквивалентной системы лежал внутри круга с радиусом, равным единице.

При выполнении достаточного условия сходимости оценка погрешности решения на k- ой итерации дается выражением:

Следует подчеркнуть, что это неравенство дает завышенное число итераций k, поэтому редко используется на практике

1.4 Численный метод решения задачи

Имеет место следующее достаточное условие сходимости метода простых итераций.

Если используется метод Якоби (выражения (1.18) для эквивалентной СЛАУ), то

достаточным условием сходимости является диагональное преобладание матрицы A, т.е.

где x^·- точное решение СЛАУ.

Процесс итераций останавливается при выполнении условия , где ее?)(kе - задаваемая вычислителем точность.

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

откуда получаем априорную оценку числа итераций k при ||б||<1

Следует подчеркнуть, что это неравенство дает завышенное число итераций k, поэтому редко используется на практике.

1.6 Текст программы

program Yakobi;

uses crt;

const

maxn=100;

type

matrix=array[1..maxn,1..maxn] of real;

vector=array[1..maxn] of real;

vector1=array[1..maxn] of real;

var

i,j,n,k1: integer;

e,norma:real;

a: matrix;

b: vector;

x2,x3: vector1;

imya,dannble_i_rezultat,ekran:string;

procedure input(var kolvo:integer; var pogreshnostb:real; var matr1:matrix; var matr2:vector);

{Ввод исходных данных}

var i,j,code:integer;

a_s: string;

b_s: string;

begin

writeln(введите количество уравнений);

readln(kolvo);

writeln(введите погрешность вычисления);

readln(pogreshnostb);

writeln(введите матрицу коэффициентов при неизвестных);

for i:=1 to kolvo do

for j:=1 to kolvo do

begin

repeat

begin

writeln(введите элемент [,i,,,j,] и нажмите Enter);

readln(a_s);

if (i=j) and (a_s=0) then

repeat

begin

writeln(Элементы на главной диагонали должны быть отличны от нуля. Повторите ввод);

readln(a_s);

end;

until a_s<>0;

val(a_s,matr1[i,j],code);

end;

until code=0;

end;

writeln(введите вектор свободных коэффициентов);

for j:=1 to kolvo do

begin

repeat

writeln(введите элемент ,j, и нажмите Enter);

readln(b_s);

val(b_s,matr2[j],code);

until code=0;

end;

1.7 Тестовый пример

2. Полиномиальная интерполяция функции методом Ньютона с разделенными разностями

Содержание