Файл: Лекция 15 Численные методы.docx

где А – матрица m m, x = (x₁, x₂,…, x_m)^T – искомый вектор, f = (f₁, f₂,…, f_m)^T – заданный вектор. Предполагается, что определитель матрицы А отличен от нуля, так что решение x существует и единственно. Для большинства вычислительных задач характерным является большой порядок матрицы А. Систему (1) можно решать по крайней мере двумя способами: или по формулам Крамера, или методом Гаусса последовательного исключения неизвестных. При больших m первый способ, основанный на вычислении определителей, требует порядка m! Арифметических действий, а метод Гаусса – только О(m³) действий. Поэтому метод Гаусса широко используется при численном решении задач линейной алгебры.

Методы численного решения системы линейных алгебраических уравнений делятся на две группы: прямые методы (Гаусса, Крамера) и итерационные методы. Итерационные методы состоят в том, что решение x системы находится как предел при nпоследовательных приближений x⁽ⁿ⁾, где n – номер итерации. Обычно задается малое число и вычисления проводятся до тех пор, пока не будет выполнена оценка

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

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

Прямые методы не предполагают, что матрица А имеет какой-либо специальный вид и ее порядок не превышает 100.

Рассмотрим систему линейных алгебраических уравнений

, (1)

где А – матрица mm, x = (x₁, x₂,…, x_m)^T – искомый вектор, f = (f₁, f₂,…, f_m)^T – заданный вектор. Предполагается, что определитель матрицы А отличен от нуля, так что решение x существует и единственно.

Запишем систему (1) в развернутом виде

(2)

…

Метод Гаусса решения системы (2) состоит в последовательном исключении неизвестных x₁, x₂,…, x_m из этой системы. Предположим, что a₁₁<>0. Поделив первое уравнение на a₁₁ получим

, (3)

где

Рассмотрим теперь оставшиеся уравнения системы (2):

. (4)

Умножим (3) на a_i₁ и вычтем полученное уравнение из i-го уравнения системы (4), i=2,…,m. В результате получим следующую систему уравнений:

(5)

Здесь обозначено

. (6)

Матрица системы (5) имеет вид

Матрицы такой структуры принято обозначать так:

Где крестиками обозначены ненулевые элементы. В системе (5) неизвестное x₁ содержится только в первом уравнении, поэтому в дальнейшем достаточно иметь дело с укороченной системой уравнений

, (7)

Тем самым мы осуществили первый шаг метода Гаусса. Если , то из системы (7) совершенно аналогично можно исключить неизвестное x₂ и прийти к системе, эквивалентной (2) и имеющей матрицу следующей структуры:

При этом первое уравнение системы (5) остается без изменения.

Исключая таким же образом неизвестные x₃,x₄,…,x_m, придем окончательно к системе уравнений вида

Матрица этой системы

(8)

. (9)

Содержит нули всюду ниже главной диагонали. Матрицы такого вида называются верхними треугольными матрицами.

Получение системы (8) составляет прямой ход метода Гаусса. Обратный ход заключается в нахождении неизвестных x₁,x₂,…,x_m из системы (8). Поскольку матрица имеет треугольный вид, можно последовательно, начиная с x_m, найти все неизвестные. Общие формулы обратного хода имеют вид

. (10)

При реализации в программе прямого хода метода Гаусса нет необходимости действовать с переменными x₁,x₂,…,x_m. Достаточно указать алгоритм, согласно которому исходная матрица А преобразуется к треугольному виду (9), и указать соответствующие преобразования правых частей системы. Получим эти общие формулы. Пусть реализованы первые k-1 шагов, т.е. уже исключены переменные x₁,x₂,…,x_k_-1. Тогда имеем систему

(11)

Рассмотрим k-е уравнение этой системы

И предположим, что . Поделив обе части этого уравнения на , получим

, (12)

где

Далее, умножим уравнение (12) на и вычтем полученное соотношение из i-го уравнения системы (11), где i=k+1,k+2,…,m. В результате последняя группа уравнений системы (11) примет вид

где

Таким образом, в прямом ходе метода Гаусса коэффициенты уравнений преобразуются по следующему правилу:

. (13)

Вычисление правых частей системы (8) осуществляются по формулам

(15)

. (16)

Коэффициенты c_ij и правые части y_i, i=1,2,..m, j=i+1,i+2,…, m, хранятся в памяти и используются при осуществлении обратного хода по формулам (10).

Основным ограничением этого является предположение о том, что все элементы , на которые проводится деление, отличны от нуля. Число называется ведущим элементом на k-м шаге исключения. Даже если какой-то ведущий элемент не равен нулю, а просто близок к нему, в процессе вычислений может происходить сильное накопление погрешностей. Выход из этой ситуации заключается в том, что в качестве ведущего элемента выбирается не

В листинге представлен метод Gauss(), который получает матрицу А и столбец свободных членов В, и реализует прямой и обратных ход метода Гаусса.

Листинг 15.1. Метод Гаусса

class Program

{

public static void Gauss(float[,] A, float[] B,

ref float[] X)

{

int Count = B.Length;

// число неизвестных = число уравнений

// решение системы: прямой ход

for (int i = 0; i <= Count - 2; i++)

// по строкам

for (int j = i + 1; j <= Count-1; j++)

// по столбцам

{

for (int k=i+1; k <= Count-1; k++)

// по строкам

A[k,j] += -A[i,j]*A[k,i]/A[i,i];

B[j] += -B[i]*A[i,j] / A[i, i];

}

// решение системы: обратный ход

for (int j = Count-1; j >= 0; j--)

// по столбцам

{

X[j] = B[j];

for (int k = j + 1; k <= Count-1; k++)

// по строкам

X[j] += -A[k, j] * X[k];

X[j] /= A[j, j];

}

Из главного модуля Main() метод Gauss() вызывается для системы

Листинг 15.2. Вызов метода Гаусса

static void Main(string[] args)

{

float [,] A = {{1,1},{1,-1}};

float[] B = { 1, 1 };

float[] X = new float[B.Length];

Gauss(A,B,ref X);

int L = X.Length;

for (int i = 0; i <= L - 1; i++)

Console.WriteLine("X[{0}] = {1}",i,X[i]);

Console.ReadKey();

}

В результате работы программы на экран выведено следующее решение:

15.2. Приближенное вычисление производных

Задача численного дифференцирования состоит в приближенном вычислении производных функции u(x) по заданным в конечном числе точек значениям этой функции. Пусть на [a,b] введена сетка:

И определены значения u_i=u(x_i) функции u(x) в точках сетки. В качестве приближенного значения u’(x_i) можно взять, например, любое из следующих разностных отношений:

В качестве примера рассмотрим функцию

для которой производная равна .

Код программы представлен в листинге 15.3.

Листинг 15.3. Приближенное вычисление производной

float a = 0;

float b = 3;

int m = 50;

float dx = (b - a) / m;

float x = 0, y = 0, yt = 0;

for (int i = 1; i <= m; i++)

{

x += dx;

y = (F3(x) - F3(x - dx)) / dx; // левая разность

yt = (F3(x + dx) - F3(x - dx)) / (2 * dx);

// центральная разность

g.DrawEllipse(Pens.Black,II(x)-2,JJ(y)-2, 4, 4);

g.DrawRectangle(Pens.Black,II(x)-2,JJ(yt)-2,4,4);

}

Результаты расчетов представлены на рис. 15.1.

Рис. 15.1. Графики функций для расчетов

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

15.3. Приближенное вычисление интегралов

Рассмотрим способы приближенного вычисления определенных интегралов

(1)

основанном на замене интеграла конечной суммой

где ck – числовые коэффициенты и xk – точки отрезка [a,b], k=0,1,…,n. Приближенное равенство

(2)

называется квадратурной формулой, а сумма вида (2) – квадратурной суммой. Точки x_k называются узлами квадратурной формулы, а числа c_k - коэффициентами квадратурной формулы.

Введем на [a,b] равномерную сетку с шагом h, т.е. множество точек

и представим интеграл (1) в виде суммы интегралов по частичным отрезкам:

(3)

Для построения формулы численного интегрирования на всем отрезке [a,b] достаточно построить квадратурную формулу для интеграла

(4)

На частичном отрезке [x_i_-1,x_i] и воспользоваться свойством (3).

Рис. 15.2. Геометрический смысл формулы прямоугольников

Формула прямоугольников. Заменим интеграл (4) выражением f(x_i_-1/2)h, где x_i_-1/2=x_i-h/2. Геометрически такая замена означает, что площадь криволинейной трапеции ABCD заменяется площадью прямоугольника ABC’D’ (рис. 6). Тогда получаем формулу

(5)

которая называется формулой прямоугольников на частичном отрезке [x_i_-1,x_i].

Рис. 15.3. Вычисление интеграла по формуле прямоугольников

Суммируя равенства (5) по I от 1 до N, получаем составную формулу прямоугольников

(6)

Формула трапеций. На частичном отрезке эта формула имеет вид

(7)

И получается путем замены подынтегральной функции f(x) интерполяционным многочленом первой степени, построенным по узлам xi-1,xi, то есть функцией

Рис. 15.4. Вычисление интеграла по формуле трапеций

Составная формула трапеций имеет вид

(8)

где .

Формула Симпсона. При аппроксимации интеграла (4) заменим функцию f(x) параболой, проходящей через точки (x_j,f(x_j)), j=i-1,i-0,5,i, т.е. представим приближенно f(x) в виде

Где - интерполяционный многочлен Лагранжа второй степени,

(9)

Проводя интегрирование, получим

Таким образом, приходим к приближенному равенству

которое называется формулой Симпсона или формулой парабол.

Рис. 15.5. Вычисление интеграла по формулам Симпсона

На всем отрезке [a,b] формула Симпсона имеет вид

Чтобы не использовать дробных индексов, можно обозначить

и записать формулу Симпсона в виде

. (10)

Все три приближенные способа представлены процедурами MethodRectangle(), MethodTrapezes() и MethodSimpson() в следующем листинге:

Листинг 15.4. Приближенное вычисление интеграла

static float F(float x)

{

return 3 * x * x;

}

static float F1(float x)

{

return x * x * x;

}

static float MethodRectangle(float a, float b)

{

int n = 200;

float dx = (b - a) / n;

float result = 0;

for (int i = 1; i <= n; i++)

result += dx * F(a+i*dx);

return result;

}

static float MethodTrapezes(float a, float b)

{

int n = 200;

float dx = (b - a) / n;

float result = (F(a)+F(b))*dx/2;

for (int i = 1+1; i <= n-1; i++)

result += dx * F(a + i * dx);

return result;

}

static float MethodSimpson(float a, float b)

{

int n = 200;

float dx = (b - a) / n;

int[] k = new int[] { 2, 4 };

float result = (F(a) + F(b)) * dx / 3;

for (int i = 1 + 1; i <= n - 1; i++)

result += k[i % 2]*dx*F(a + i * dx)/3;

return result;

}

static void Main(string[] args)

{

float a = 0;

float b = 1;

Console.WriteLine("Точное={0}",F1(b)-F1(a));

float s = MethodRectangle(a, b);

Console.WriteLine("Метод прямоугольников");

Console.WriteLine("s = {0}", s);

Console.WriteLine("");

s = MethodTrapezes(a, b);

Console.WriteLine("Метод трапеций");

Console.WriteLine("s = {0}", s);

Console.WriteLine("");

s = MethodSimpson(a, b);

Console.WriteLine("Метод Симпсона");

Console.WriteLine("s = {0}", s);

Console.WriteLine("");

Console.ReadKey();

}

В качестве примера рассмотрим определенный интеграл:

Полученные результаты показывают следующее:

15.4. Линейные дифференциальные уравнения

Рассмотрим задачу Коши для системы обыкновенных дифференциальных уравнений

(1)

или подробнее

(2)

(3)

Существует две группы численных методов решения задачи Коши: многошаговые разностные методы и методы Рунге-Кутта. Приведем примеры и поясним основные понятия, возникающие при использовании численных методов. Для простоты будем рассматривать одно уравнение

. (4)

Введем по переменному t равномерную сетку с шагом >0, т.е рассмотрим множество точек

Будем обозначать через u(t) точное решение задачи (4), а через y_n=y(t_n) – приближенное решение. Заметим, что приближенное решение является сеточной функцией, т.е. определено только в точках сетки .