Основы вычислительной математики

Ю.П. Хрусталев

 

 

ОСНОВЫ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ

Практикум

Для студентов, обучающихся по специальностям:

: 2201 “Вычислительные машины, комплексы, системы и сети ”

2202 “Автоматизированные системы обработки информации и управления ”

719 «Информационные системы в технике и технологиях»

 

 

Издательство

Иркутского государственного технического университета

Рецензент – начальник отдела ИСЭМ СО РАН, д.т.н. профессор В.И. Зоркальцев

Основы вычислительной математики. Практикум. Разработал доцент кафедры вычислительной техники Ю.П. Хрусталев – Иркутск, 2008 – с 94

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

Приведены алгоритмы решения основных задач вычислительной математики и показана методика их реализации с помощью системы MATHCAD.

Для студентов, обучающихся по специальностям:

: 2201 “Вычислительные машины, комплексы, системы и сети ”

2202 “Автоматизированные системы обработки информации и управления ”

719 «Информационные системы в технике и технологиях»

 

Библиография – Назв. – 7 Табл. – 6 http://detc.usu.ru/Assets/aMATH0021/l1_3.htm

ОГЛАВЛЕНИЕ

 

1. Раздел 1. Погрешности вычислений Стр. 4

2. Раздел 2. Численные методы решения систем линейных

алгебраических уравнений Стр. 14

3. Раздел 3. Решение нелинейных уравненийСтр. 40

 

4. Раздел 4. Интерполирование функций Стр. 46

 

5. Раздел 5. Численное интегрирование Стр. 65

 

6. Раздел 6. Численные методы решения обыкновенных

дифференциальных уравнений Стр. 85

Раздел 1

Погрешности вычислений

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

Абсолютная погрешность ε - это разность между истинным и приближенным значением величины:

(1)

Относительная погрешность - это отношение величины абсолютной погрешности к приближению:

(2)

Для величин, близких по значению к единице, абсолютная и относительная погрешности имеют один порядок. Для очень больших и очень маленьких величин эти погрешности существенно отличаются. Например, если точное значение некоторой величины равно 0.00006, а приближенное значение 0.00005, то = , а ( 20%). Если точное значение величины равно 100500, а приближенное 100000, то , а .

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

 

СЛОЖЕНИЕ (ВЫЧИТАНИЕ)

(3)

УМНОЖЕНИЕ

(4)

ДЕЛЕНИЕ

(5)

Для уменьшения погрешности в практической деятельности рекомендуется соблюдать следующие правила:

1) По необходимости избегать вычитания двух почти равных чисел.

2) Использовать формулы типа

a(b-c) = ab -ac; (b-c)/a = b/a - c/a.

 

Причины возникновения погрешностей

1) Неточное представление исходной информации.

2) Округление.

3) Ограничение бесконечного математического процесса конечным числом операций.

Задача лабораторной работы - определение погрешностей третьего вида.

 

 

Лабораторная работа 1

Определение погрешностей при вычислении функций методом разложения их в степенной ряд

Работа выполняется с использованием палитры программирования системы автоматизации математических вычислений Mathcad и системы Excel.

 

Задание на работу: вычислить значение функции по приближенной формуле с заданной точностью:

1) получить рекуррентное соотношение для вычисления последующего члена степенного ряда по предыдущему;

2) вычислить значение функции с помощью программы на Mathcadе;

3) вычислить значение функции, используя систему EXEL.

 

Варианты заданий

	Функция f (x)	Сумма (приближенное значение) Q (x)	Значение аргумента x	Точ-ность ()
	Cos(x)		0.5
	(x)		0.8
	(x)		0.5
			0.5
			3.5
			1.5
			0.7
			0.5
			0.5
	=2Arcth(x)		0.5
	=2Arcth(x)		2.5
			3.5
			3.5
			2.5
			0.6
			0.5
			0.5
			0.5
			0.8
			0.7
			0.8

 

Требования к оформлению отчета

Отчет по лабораторной работе должен содержать:

- название работы

- цель работы

- математическую формулировку задачи (представление заданной функции степенным рядом, вывод рекуррентного соотношения и алгоритм решения)

- тексты программы (в системе Mathcad и Excel)

- результаты, полученные в процессе выполнения работы

- выводы

Вспомогательные материалы

Пример вычисления функции f(x) = Sin(x) с заданной точностью

eps = 0.000001

X = 0.5

1. Алгоритм решения задачи

 

1.1 В сумму, предварительно обнуленную, заносится первый член ряда.

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

1.3 С помощью рекуррентной формулы вычисляется следующий член ряда и добавляется к сумме.

1.4 Переход на пункт 1.2.

2. Вывод рекуррентного соотношения

(При выводе данного соотношения в общем члене ряда индекс i заменяется на i+1, полученный i+1 -й член делится на предыдущий)

 

3. Программа для вычисления синуса

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

 

Program Lsb1;

{****************************************************************

ПОГРЕШНОСТИ ВЫЧИСЛЕНИЙ

****************************************************************}

Uses Crt;

Var

x,e,eps,a0,a1,r,s,y,z: real;

i,n,k: integer;

BEGIN

ClrScr;

WRITE(' Введите Х = ');

READLN (x);

WRITE (' Введите е = ');

READLN (e);

y:= Sin(x);

WRITELN (' SIN (x) = ', y);

WRITELN;

a0:= x;

i:= 0;

eps:= y - a0;

WRITELN (' N', ' S ', ' E');

WRITELN (' ', i+1, ' ', a0, ' ', eps);

r:= - x*x;

i:= 1;

s:= a0;

REPEAT

n:= 2*i;

z:= n*(n+1);

a1:= r *a0/z;

s:= s+a1;

eps:= y-s;

WRITELN (' ', i+1, ' ', s, ' ', eps);

a0:= a1;

i:= i+1;

UNTIL abs (eps) < e;

END.

Результаты работы программы

Введите X = 0.5

Введите е = 0.0000001

SIN (x) = 4.7942553860E-01

 

N S E

 

1 5.0000000000E-01 -2.0574461396E-02

2 4.7916666667E-01 2.5887193760E-04

3 4.7942708333E-0.1 -1.5447290025E -06

4 4.7942553323E-0.1 5.3701114666E- 09

 

4. Пример вычисления синуса в системе Excel

 

	пример вычисления SIN(X) в Exсel
	текущий член суммы	сумма	точное знач.	погрешность	x * x
	"-x^2/((2*i)(2*i+1)) * A[i-1]	S[i]= S[i-1] + A[i]	Sin(x)	Sin(x) - S(i)	x^2
	Формулы, занесенные в соответствующие столбцы второй строки (строка 1)
	"=-($F$11)*B11/(2 *A12*(2*A12+1))"	"=C11+ B12"		"=$D$11 - C12"
	0,5	0,5	0,479425539	-0,020574461	0,25
	-0,020833333	0,479166667		2,58872E-04
	2,60417E-04	0,479427083		-1,54473E-06
	-1,55010E-06	0,479425533		5,37008E-09
	5,38229E-09	0,479425539		-1,22129E-11
	-1,22325E-11	0,479425539		1,95954E-14
	1,96033E-14	0,479425539
	-2,33373E-17	0,479425539

5. Решение задачи в системе Mathcad с использованием палитры программирования

5.1 Разработка программы в MathCAD

Чтобы раскрыть палитру программирования, следует выполнить действия, показанные на рисунке 1.1

Рис. 1.1 Обращение к палитре программирования

 

Оператор Add Line обозначает в MathCAD программный модуль - вертикальную черту, справа от которой последовательно записываются операторы.

Оператор ç означает присваивание: F ç 2+B

Otherwise – это оператор аналогичный Else в Pascal. С остальными операторами вы знакомы.

Создадим новый лист в MathCAD. Напишем оператор ORIGIN:= 1 – для того, чтобы MathCAD начинал нумерацию в массивах с 1, а не с нуля (по умолчанию нумерация начинается с нуля). Затем вводим исходные данные

Рис 3.

Напишем имя функции, оператор присваивания и вставим программный модуль, несколько раз нажав кнопку Add Line на панели программирования. Затем начинаем набирать программу. Пример программы приведен на рис. 1.2

 

Рис 1.2 Пример программы для вычисления функции y = Sin(x)

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

Оператор while содержит сложное условие выхода из цикла. Помимо проверки достижения заданной точности ограничено также число итераций (в приведенном примере число итераций не должно превышать пяти, хотя в общем случае оно должно быть много больше). Это ограничение позволяет избежать «зацикливания» программы, если она не верна, и процесс приближения функции степенным рядом не сходится.

Текущие значения суммы и членов ряда запоминаются в одномерных массивах SS и B. Это необходимо для отладки программы. Кроме того, значения текущей суммы и членов степенного ряда позволяют наглядно представить процесс приближения функции.

Программа возвращает переменную, расположенную в последней строке программы. Если необходимо в результате работы программы получить значения нескольких переменных, то их следует поместить в скобки, т. е. вывести вектор. Вектор задается с помощью палитры Math (изображение матрицы). При этом число строк должно равняться числу возвращаемых переменных. Матрица содержит один столбец.

На рис.1.3 приведен пример вывода на экран результатов работы программы. В верхней части рисунка показан полный вектор возвращаемых переменных. Третья и четвертые строки сами являются векторами, на что указывает их изображение (фигурные скобки с указанием размерности выводимых матриц – 2,1, т. е. вектор, содержащий две строки). В данном примере число итераций равнялось двум. Ниже показан способ отображения этих векторов.

Рис. 1.3 Вывод на экран дисплея результатов работы программы

 

Контрольные вопросы

1.Назовите причины возникновения погрешностей.

2.Может ли относительная погрешность превышать абсолютную погрешность?

3.В каком случае относительная погрешность меньше абсолютной погрешности?

4.В каких единицах выражаются абсолютная и относительная погрешности?

 

Раздел 2

Варианты заданий.

( Система уравненийв матричной записи имеет вид Ax = b)

1. A= , b= 2. A= , b=

 

3. A= , b= 4. A= , b=

 

5. A= , b= 6. A= , b=

7. A= , b= 8. A= , b=

9. A= , b= 10. A= , b=

 

11. A= , b= 12. A= , b=

13. A= , b = 14. A= , b=

 

15. A= , b = 16. A= , b=

 

17. A= , b = 18. A= , b=

 

19. A= , b = 20. A= , b=

 

21. A= , b =

 

Требования к оформлению отчета

Отчет должен содержать:

- Название и цель работы

- Задание на работу

- Текст программы на Mathcadе

- Результаты работы программы

- Проверку решения

 

Вспомогательные материалы

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

Этап 1. Объявление переменных

Определим переменные и зададим матрицу системы линейных уравнений:

Создадим новый лист в MathCAD и для начала напишем оператор ORIGIN:= 1 – это для того чтобы MathCAD начинал нумерацию в матрицах с 1, а не с нуля.

Теперь определим матрицу А.

Begin

d1:=0;

da:=0;

lin:=False;

For i:=0 to n-2 do

For j:=i+1 to n-1 do

Begin

For k:=0 to n-1 do

If matr[j,k]<>0 then d1:=d1+matr[i,k]/matr[j,k];

d1:=d1/n;

For k:=0 to n-1 do

If matr[j,k]<>0 then If d1=matr[i,k]/matr[j,k] then da:=da+1;

If da=n then lin:=True;

d1:=0;

da:=0;

end;

If lin then Proverka:=True Else Proverka:=False;

end;

Лабораторная работа 3

Раздел 3

Метод половинного деления

Описание метода.

Пусть дано уравнение

f(x) = 0, (1)

причем функция f(x) непрерывна на отрезке [ a, b] и f(a)f(b) < 0.

Для вычисления корня уравнения (1), принадлежащего отрезку [ a, b], найдем середину этого отрезка x₁ = (a + b) / 2. Если f(x₁) 0, то для продолжения вычислений выберем ту из частей банного отрезка [ a, x₁] или [ x_1, b], на концах которой функция f(x) имеет противоположные знаки. Концы нового отрезка обозначим через a₁ и b_{1 .}

Новый суженный промежуток [ a₁, b₁ ] снова делим пополам и продолжаем вычисления по разработанной схеме и т. д. В результате получаем либо точный корень уравнения (1) на каком - то этапе, либо последовательность вложенных отрезков [ a, b ], [ a₁, b₁ ], ..., [ a _n, b _n ],.. таких, что

f(a _n)f(b _n) < 0 (n = 1, 2,...), (2)

b _n - a _n = (1/ 2 ⁿ) (b- a) (3)

Число ξ - общий предел последовательности { a _n } и { b _n } - является корнем уравнения

f(x) = 0.

Оценку погрешности решения на n -ом шаге вычислений можно получить из соотношения (3) в виде

 

0 < ξ - a _n (1/ 2 ⁿ) (b- a) = b _n - a _n (4)

Здесь a _n ξ c точностью ε не превышающей (1/ 2 ⁿ) (b- a).

Y

		Y = f (x)

b

X

 

 

 

Рис. 3.1 Наличие единственного корня уравнения на интервале [a,b]

 

Блок – схема алгоритма, реализующего метод половинного деления, приведена на рис. 3.3.

Теорема.

Пусть функция g (x) имеет на отрезке [a, b] непрерывную производную и выполнены два условия:

1) q < 1 при x [a, b];

2) значения функции y = g(х) принадлежат отрезку [a,b] для любого x [a, b]

Тогда при любом выборе начального приближения x⁽⁰⁾ [a, b] процесс итераций сходится к единственному корню уравнения (1) на отрезке [a, b]

Оценка погрешности k -го приближения x ^(k) к корню такова:

, (8)

где

Укажем теперь один из способов преобразования уравнения

f(x) = 0 (9)

к виду x = g(x), допускающему применение метода итераций, сходящихся к решению уравнения (9).

Для любого числа уравнение (9) равносильно уравнению (5), где

g (x) = x + f(x).

Предположим, что производная f ^' (x) > 0 и непрерывна

на [ a,b ]. Пусть , ;

положим

,

и рассмотрим функцию

. (10)

Для функции, определенной формулой (10), выполняются достаточные условия сходимости метода итераций решения уравнения (9). В частности, условие 1) теоремы следует из неравенств

0 < m f ^' (x) M,

0 g ^' (x) = 1 - (1/M) f ^' (x) 1 - m/M = g < 1 .

Замечание1. Если окажется, что производная f ^' (x) отрицательна на отрезке [ a, b], то уравнение (1) можно заменить на равносильное уравнение -f(x) = 0 и использовать указанное преобразование.

Замечание 2. Если вычисление точного числа затруднительно, то можно заменить его произвольным числом М₁> M. Однако при большом М₁ число q = 1 - m / М₁ ближе к единице и процесс итераций сходится медленнее.

Замечание 3. При нахождении корня уравнения (1) с заданной точностью или при оценке погрешности k-го приближения можно, не вычисляя точного значения числа

q = max | g ^' (x) |,ограничиться следующей практической рекомендацией:

при 0 < q (1/2) (11)

при (1/2) < q < 1. (12)

Блок – схема алгоритма, реализующего итерационный метод, приведена на рис. 3.2.

 

 

>

<

<

b = x1

d = b - a

a = x1

0

x1 = (a + b ) / 2

 

 

 

Рис 3.3 Блок – схема алгоритма, реализующего метод половинного деления

 

Лабораторная работа 4

Варианты заданий.

1. x ⁴ - 3x -20 = 0 (x > 0) 2. x ³ - 2x - 5 = 0 (x > 0)

3. x ³ + 3x + 5 = 0 4. x ⁴ + 5x -7 = 0 (x > 0)

5. x ³ - 12x - 5 = 0 (x > 0) 6. x ³ - 2x ² - 4x + 5 = 0 (x < 0)

7. x + e ^x = 0 8. x ⁵ - x - 2 = 0

9. x ³ - 10x + 5 = 0 (x < 0) 10. 2 - lnx - x = 0

11. x ³ + 2x - 7 = 0 12. x ³ + x ² - 11 = 0 (x > 0)

13..x ⁴ - 2x - 4 = 0 (x > 0) 14. 2e ^x + x - 1 = 0

15. x ⁴ - 2x - 4 = 0 (x < 0) 16. 2x ³ + x ² - 4 = 0 (x > 0)

17. e ^x - x - 2 = 0 18. (1/2) e ^x - x - 1 = 0 (x > 0)

19. x ² - cos x = 0 (x > 0) 20. x ² + lnx = 0

Требования к оформлению отчета

Отчет по лабораторной работе должен содержать:

- название работы

- цель работы

- тексты программ

- результаты, полученные в процессе выполнения работы

- выводы

Вспомогательные материалы

Для разработки программ на Mathcad е можно использовать приемы, описанные в лабораторных работах №1 и №2.

Существует бесчисленное множество линейных функций, имеющих корень на интервале[a,b]. Поясним это на примере.

Пусть a = 2, b = 7. Пусть корень уравнения равен 5. Тогда функция

F(x) = k1 + k2*x, имеющая корень, равный пяти, может иметь такой вид

F(x) = 10 – 2*x (один из коэффициентов задается произвольно, другой находится из уравнения F(x) = 0).

 

Контрольные вопросы:

1. Зачем нужна процедура отделения корней?

2. Что называется корнем уравнения?

3. Какова точность метода половинного деления?

4. Каким образом исходное уравнение преобразуется к виду, удобному для итераций?

5. Чему равна оценка погрешности k -го приближения?

Раздел 4

Интерполирование функций

Краткое введение. Постановка задачи интерполирования.

На отрезке заданы n+1 точки , которые называются узлами интерполяции, и значение некоторой функции в этих точках

. (1)

Требуется построить интерполирующую функцию F(x), принадлежащую известному классу и принимающую в узлах интерполяции те же значения, что и , т.е.

(2)

 

Рис. 4.1 Интерполирование функции y = f(x)

 

В общем случае, задача имеет бесчисленное множество решений.

Задача становится однозначной, если решение искать в заданном классе функций.

Будем искать полином степени не выше n и удовлетворяющий условию (2).

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

Если - имеет место задача интерполирования (интерполирование “в узком смысле”).

При решается задача экстраполирования.

 

Обратное интерполирование

Задача обратного интерполирования: по заданному значению функции найти аргумент , при котором . Функция y=f(x) задана таблично.

Предположим, что на отрезке [a, b], содержащем узлы интерполяции, функция f(x) монотонна. Тогда существует однозначная обратная функция x=F(y). Она задана той же таблицей, что и y=f(x), только теперь аргументом будет значение , а -соответствующее значение функции.

В этом случае обратное интерполирование сводится к обычному интерполированию для функции x=F(y). Т.е. строится интерполяционный многочлен (например, по формуле Лагранжа) – многочлен . При подстановке в значения - получаем .

Второй способ применим ко всякой функции f(x) (не обязательно к монотонной). Не меняя ролями функцию и аргумент, записываем по какой – либо формуле интерполяционный многочлен