Вычислительная математика (стр. 6 )

(define (int fun points coefs h)

(* h (apply + (map * coefs (map fun points)))))

; Подынтегральная функция

(define (fx x) (* (sqrt x) (cos x)))

; Формирование списка узлов

(define (make-point-list a h n)

(cond

((= 0 n) '())

(else (cons a (make-point-list (+ a h) h (- n 1))))))

; Списки весов для левых/правых/средних прямоугольников

(define (mid n) (if (= 0 n) '(0) (cons 1 (mid(- n 1)))))

(define (left n) (if (= 0 n) '(0) (cons 1 (left (- n 1)))))

(define (right n) (cons 0 (mid (- n 1))))

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

(define (parabol n) (cons (/ 1 3) (parabol1 (- n 1))))

(define (parabol1 n)

(cons (/ 4 3) (parabol2 (- n 1))))

(define (parabol2 n)

(if (<= n 1)

(list (/ 1 3))

(cons (/ 2 3) (parabol1 (- n 1)))))

(define (trap n) (cons (/ 1 2) (trap1 (- n 1))))

(define (trap1 n) (if (= 0 n) '(1/2) (cons 1 (trap1 (- n 1)))))

; Подсчёт списка значений интеграла всеми методами с заданным числом узлов

(define (iN n)

(let ((h (/ 1 n)) (points (make-point-list 0 (/ 1 n) (+ n 1))))

(list (int fx points (left n) h)

(int fx points (right n) h)

(int fx (map (lambda (x) (+ x (/ h 2))) points) (mid n) h)

(int fx points (trap n) h)

(int fx points (parabol n) h))))

; Вывод результатов

(print (iN 12))

(print (iN 24))

(print (map / (map abs (map - (iN 12) (iN 24))) '(3 3 3 3 15)))

(print (map + (iN 24)

(map / (map abs (map - (iN 12) (iN 24))) '(3 3 3 3 15))))

5.5 Вывод

Искомый определённый интеграл получился равным

– по формуле трапеций

– по формуле прямоугольников

– по формуле Симпсона

Как видно, наименее точной из этих формул оказалась формула прямоугольников, наиболее точной – формула Симпсона.

6 Задание №7

6.1 Задание

Вычислить интеграл по формуле Гаусса и оценить погрешность для :

6.2 Описание метода решения задачи

Воспользуемся следующей формулой: .

Где – нули полинома Лежандра четвёртой степени .

Коэффициенты и точки – табличные величины.

6.3 Оценка погрешности

Остаточный член формулы Гаусса для выражается следующим образом:

[Далее следует привести нахождение десятой производной и её максимума.]

6.4 Код программы на языке Visual Basic

Module Module1

Private Function f(t As Double) As Double

f = 1 / (10 + t)

End Function

Sub Main()

Dim integral As Double = 0.0

Dim t()={System. Math. Sqrt((3-2*System. Math. Sqrt(1.2))/7), _

- System. Math. Sqrt((3-2*System. Math. Sqrt(1.2))/7), _

System. Math. Sqrt((3+2*System. Math. Sqrt(1.2))/7), _

-System. Math. Sqrt((3+2*System. Math. Sqrt(1.2))/7)}

Dim A()={(18+System. Math. Sqrt(30))/36, _

(18+System. Math. Sqrt(30))/36, _

(18-System. Math. Sqrt(30))/36, _

(18-System. Math. Sqrt(30))/36}

For i = 0 To 3 Step 1

integral = integral + A(i) * f(t(i))

Next

Console. WriteLine(integral)

End Sub

End Module

6.5 Вывод

Искомый интеграл получился равным

Данный метод оказался более точным, чем другие рассмотренные методы численного интегрирования.

7 Задание №12

7.1 Задание

Сравнить эффективность метода Ньютона и метода обратной квадратичной интерполяции для достижения заданной точности. .

на промежутке .

7.2 Описание метода решения задачи

Метод Ньютона поиска корней уравнения был описан выше.

Метод обратной квадратичной интерполяции заключается в следующем:

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

Продолжаем этот процесс, пока .

Таблица 8 – Сравнение эффективности метода Ньютона и метода обратной квадратичной интерполяции

Метод Ньютона	Метод обратной квадратичной интерполяции
0	0,2000	0	0,4400
1	0,5423	1	0,4384
2	0,4407	2	0,4384
3	0,4384

7.3 Код программы на языке C++

#include <iostream>

#include <fstream>

#define _USE_MATH_DEFINES

#include <cmath>

float f(float q) {

return q * q - cos(M_PI * q);

}

float f1(float q) {

return 2 * q + M_PI * sin(M_PI * q);

}

int main () {

std::ofstream out("z12.out");

std::out. precision(4);

std::out. setf (ios::fixed);

std::out << "Square:" << std::endl;

float x0 = 0.2f, x1 = 0.4f, x2 = 0.5f;

float y0 = f(x0), y1 = f(x1), y2 = f(x2);

float e = 0.0001f;

float x = (y0 * y0 * y1 * x2 + y2 * y2 * y0 * x1 + y1 * y1 * y2 * x0

- y2 * y2 * y1 * x0 - y0 * y0 * y2 * x1 - y1 * y1 * y0 * x2)

/ (y0 * y0 * y1 + y2 * y2 * y0 + y1 * y1 * y2 - y1 * y2 * y2

- y2 * y0 * y0 - y0 * y1 * y1);

int i = 0;

while (fabs(x - x2) >= e)

{

std::out << i << ' ' << x0 << ' ' << x1 << ' ' << x2

<< ' ' << x << std::endl;

x0 = x1; x1 = x2; x2 = x;

y0 = f(x0); y1 = f(x1); y2 = f(x2);

x = (y0 * y0 * y1 * x2 + y2 * y2 * y0 * x1 + y1 * y1 * y2 * x0

- y2 * y2 * y1 * x0 - y0 * y0 * y2 * x1 - y1 * y1 * y0 * x2)

/ (y0 * y0 * y1 + y2 * y2 * y0 + y1 * y1 * y2 - y1 * y2 * y2

- y2 * y0 * y0 - y0 * y1 * y1);

i++;

}

std::out << i << ' ' << x0 << ' ' << x1 << ' ' << x2

<< ' ' << x << std::endl;

std::out << x << std::endl;

std::out << "Newton:" << std::endl;

x0 = x1 = 0.2f;

x1 = x0 - f(x0) / f1(x0);

i = 0;

while (fabs(x0 - x1) >= e)

{

out << i << ' ' << x0 << ' ' << x1 << endl;

x0 = x1; x1 = x0 - f(x0) / f1(x0);

i++;

}

std::out << i << ' ' << x0 << ' ' << x1 << std::endl << x1;

return 0;

}

7.4 Вывод

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

8 Задание №13

8.1 Задание

Вычислить коэффициенты обусловленности корней полинома пятой степени в зависимости от коэффициентов полинома a priori и a posteori, то есть теоретически и вычислив корни полинома с помощью приближенных методов по формуле Ньютона.

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

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

8.2 Алгоритм определения машинной точности

printf( "\tЭпсилон\tРазмер мантиссы\n");

float epf = 1;

int i = 1;

while ((1+epf/2) != 1)

{

epf/=2;

++i;

}

printf( "float:\t%G\t%i\n", epf/2, i);

Величина машинного эпсилон для типа float составляет , длина мантиссы – 24.

8.3 Априорная оценка погрешности

Вычислим корни уравнения.

– корень первой кратности.

– корень второй кратности.

– корень второй кратности.

Сначала найдём коэффициенты обусловленности для корня первой кратности.

Найдём коэффициенты обусловленности для корня [необходимо привести явную формулу для вычисления коэффициентов]: