">
Математика Математический анализ
Информация о работе

Тема: Вычислительная математика

Описание: Постановка. Сведение краевой задачи к задаче Коши. Решение для дифференциального уравнения. Метод достижения заданной точности. Описание результатов. Конечные разности. Решение системы уравнений методом прогонки Достижение заданной точности.
Предмет: Математика.
Дисциплина: Математический анализ.
Тип: Курсовая работа
Дата: 31.08.2012 г.
Язык: Русский
Скачиваний: 20
Поднять уникальность

Похожие работы:

Курсовая работа

По курсу: «Вычислительная математика»

Вариант №47

2004г.

Оглавление:

Постановка задачи ……..……………………………………….3

Сведение краевой задачи к задаче Коши……………………...4

Сведение краевой задачи к задаче Коши……………………………………..4

Решение задачи Коши для дифференциального уравнения…………………4

Метод достижения заданной точности……………………………………….5

Описание решения и программы……………………………………………...5

Описание результатов………………………………………………………….6

Листинг программы……………………………………………………………13

Метод конечных разностей…………………………………….7

Метод конечных разностей……………………………………………………7

Решение системы уравнкний методом прогонки…………………………….8

Метод достижения заданной точности……………………………………….8

Описание решения и программы………………………………………………8

Описание результатов…………………………………………………………..9

Листинг программы……………………………………………………………14

Метод Галёркина……………………………………………….10

Метод Галёркина………………………………………………………………10

Решение системы уравнений………………………………………………….11

Метод достижения заданной точности……………………………………….11

Описание решения и программы……………………………………………...11

Описание результатов………………………………………………………….12

Листинг программы……………………………………………………………16

Выводы…………………………………………………………..18

Литература………………………………………………………18

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

Линейным дифференциальным уравнением второго порядка с заданными краевыми условиями называется уравнение вида:

 (1.1)

(1.2)

где функции P(x), Q(x), F(x) непрерывны, и - заданные постоянные, причем:  и .

В задаче требуется найти решение y(x) уравнения (1.1), удовлетворяющее краевым условиям (1.2).

В данном варианте курсовой работы необходимо:

Решить краевую задачу для обыкновенного дифференциального уравнения

y-xy-3,1xy=3x+2

y(1,6)=-2; y(2,1)=-3; с точностью E= 1.e-5.

Проверить достигнутую точность и результаты представить с шагом h=0,05.

Решение провести тремя следующими методами:

1). Сведением краевой задачи к задаче Коши;

2). Методом конечных разностей;

3). Методом Галеркина.

Из условий следует, что: функции, P(x)=-x; Q(x)=-3,1x; F(x)=3x+2 и постоянные.

(1=1; (2=0; (1=1; (2=0; a=1,6; b=2,1; A=-2; B=-3.

II.Сведение краевой задачи к задаче Коши

2.1. Сведение краевой задачи к задаче Коши

Будем искать решение уравнения (1.1) с заданными краевыми условиями (1.2) в виде: (2.1)

Функция y(x) является линейной комбинацией функций U(x) и V(x), а C=const.

Подставив (2.1) в дифференциальное уравнение (1.1) получим систему из двух уравнений:

CU+V+P(x)[CU+V]+Q(x)[CU+V]=F(x)

C(U+P(x)U+Q(x)U)+V+P(x)V+Q(x)V=F(x)

Требуем, чтобы скобка (U+P(x)U+Q(x)U) равнялася 0

 (2.2)

где U(x)- ненулевое решение соответствующего однородного уравнения, а V(x) - некоторое решение данного неоднородного уравнения.

Первое краевое условие (1.2) должно выполняться для функции y при любом C. Используя это, получим следующее:

 (2.3)

Для того, чтобы равенство (2.3) было справедливо при любом С, необходимо и достаточно:

 (2.4)

Используя второе краевое условие (1.2) вычислим постоянную С:

,

откуда: , где 

В результате задача с краевыми условиями сводится к решению задачи Коши для функций U(x),V(x). И исходя из условий курсовой работы () получим:

,где ,  (2.5),(2.6)

и 

2.2.Решение задачи Коши для дифференциального уравнения

Рассмотрим процесс нахождения решения для уравнения (2.6). Решение уравнения (2.5) производится аналогичным способом, если в решении уравнения (2.6) учесть, что .

Согласно выбранному методу сделаем замену 

получим систему уравнений:



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

, где h - заданный шаг.

Значения  и  для функций V и z определяются формулами:

 ,где:



2.3. Метод достижения заданной точности

Для достижения заданной точности используется следующий метод:

Вычисляем y(x) два раза: один раз с шагом h, другой с шагом меньшим h (в приведенных вычислениях h/2),получая при этом значения более точные. Сравниваем значения в первом и во втором случаях. Если максимальная погрешность (норма разности между двумя полученными соответствующими значениями) не превышает заданной точности E, то выбранный шаг h можно считать достаточным и полученная функция y(x) удовлетворяет заданной точности и является ответом для данной задачи. В противном случае продолжают уменьшать шаг h, пока не будет достигнута заданная точность.

2.4. Описание результатов

Программа выдает следующие результаты:

max_eps: = 4,7e-06

x=1.60 y(x)=-2.000000 eps=0.00e+00

x=1.65 y(x)=-2.018973 eps=6.45e-07

x=1.70 y(x)=-2.048391 eps=1.26e-06

x=1.75 y(x)=-2.090075 eps=1.84e-06

x=1.80 y(x)=-2.146270 eps=2.35e-06

x=1.85 y(x)=-2.219749 eps=2.75e-06

x=1.90 y(x)=-2.313942 eps=2.99e-06

x=1.95 y(x)=-2.433102 eps=2.97e-06

x=2.00 y(x)=-2.582512 eps=2.59e-06

x=2.05 y(x)=-2.768745 eps=1.68e-06

x=2.10 y(x)=-3.000000 eps=0.00e+00

Решено с шагом h=0.025 и максимальной точностью = 2.99e-06

В полученной таблице приведено решение заданного дифференциального уравнения 2-ого порядка. В столбце х приведено разбиение отрезка [1,6;2,1] с шагом h=0.05. В столбце y соответствующие значения функции y(xn) (n=0,1,2,…,15). Погрешность вычисленных значений функции определяется максимальным значением столбца eps: max_eps, которое должно быть меньше заданной точности E= 0.00001.

III. Метод конечных разностей

3.1. Метод конечных разностей

Данный метод решения предполагает сведение краевой задачи к системе конечно-разностных уравнений. Для этого нужно разбить отрезок [a,b] на n равных частей длины h , . Точки разбиения имеют абсциссы:  (i=0,1,2,…,n),  .

В формулах приведенных ниже введены следующие обозначения: , , ,   

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

(i=1,2,3,…,n-1),(3.1)

где  - погрешность формулы порядка .

Подставив в дифференциальное уравнение (1.1) выражения (3.1), для точек  (i=1,2,3,…,n-1) его приближенно можно заменить линейной системой уравнений:

, (i=1,2,...,n-1) (3.2)

Или приведя подобные члены:

, где (i=1,2,...,n-1),(3.3)

где  ;  (3.4)

Используя краевые условия (1.2), и введя фиктивные точки  и , получим:



И исходя из условий курсовой работы () получим:

, где (i=1,2,...,n-1),(3.5)

3.2. Решение системы уравнений методом прогонки

Система уравнений (3.5) представляет линейную систему из n+1 уравнений с n+1 неизвестными , которые являются значениями искомой функции y(x) в точках . Наиболее простым методом решения таких систем является метод прогонки, имеющий порядок сложности n.

Решая уравнение (3.3) относительно уi, находим:

уi= fi/mi- ni(уi-1/mi- уi+1/mi (3.6)

предположим, что с помощью полной системы (3.3) из уравнения исключен член, содержащий уi-1. Тогда уравнение (3.6) может быть записано в виде:

уi = ci(di - уi+1) (3.7)

где ci, di – некоторые коэффициенты (i=1,2,3,…..,n-1). Найдем формулы для этих коффициентов.

Выразим из формулы (3.7) уi-1 получим:

уi-1 = ci-1(di-1 - ci-1(уi (3.8)

подставляя это выражение в формулу (3.6), будем иметь:

уi= fi/mi - ni(( ci-1(di-1 - ci-1(уi)/mi - уi+1/mi (3.9)

разрешая полученное уравнение относительно уi, находим:

уi = (fi - ni(ci-1(di-1 - уi+1)/(mi - ni(ci-1) (3.10)

Отсюда, сравнивая формулы (3.7) и (3.10), получаем:

 (i=1,2,3,...,n) (3.11)

Из первого уравнения системы (3.5) получаем.

Формулы (3.11) представляют собой прямой ход метода прогонки, в результате которого будут получены значения .

Теперь, используя формулу (3.7) и первое краевое условие, мы можем последоваельно найти уi-1, уi-2,……, у0 (обратный ход)

3.3. Метод достижения заданной точности

Для достижения заданной точности используется метод аналогичный методу рассмотренному в пункте 2.3. Точность достигается за счет уменьшения шага.

3.4. Описание результатов

Программа выдает следующие результаты:

Решение краевой задачи для дифференциального уравнения

II-го порядка с точностью Е = 0.00001 методом конечных разностей

x=1.60 y(x)=-2.000000 eps=0.00e+00

x=1.65 y(x)=-2.018975 eps=1.32e-06

x=1.70 y(x)=-2.048393 eps=2.48e-06

x=1.75 y(x)=-2.090078 eps=3.47e-06

x=1.80 y(x)=-2.146273 eps=4.22e-06

x=1.85 y(x)=-2.219752 eps=4.70e-06

x=1.90 y(x)=-2.313945 eps=4.84e-06

x=1.95 y(x)=-2.433105 eps=4.55e-06

x=2.00 y(x)=-2.582514 eps=3.73e-06

x=2.05 y(x)=-2.768746 eps=2.26e-06

x=2.10 y(x)=-3.000000 eps=0.00e+00

Результаты представлены с шагом h= 0.0012500

максимальная точность Е=4.84e-06

В полученной таблице приведено решение заданного дифференциального уравнения 2-ого порядка. В столбце х приведено разбиение отрезка [3.0;3.6] с шагом h=0.04. Погрешность вычисленных значений функции определяется максимальным значением столбца eps: max_eps, которое должно быть меньше заданной точности E= 0.0001. Было получено max_eps = 4.3e-05, что меньше заданной точности E= 1.e-4. Следовательно, при шаге h=H/4=0.01 получились результаты, удовлетворяющие условиям.

IV. Метод Галеркина

4.1. Метод Галеркина

В данном методе решения вводятся следующие обозначения:

линейный дифференциальный оператор ,

линейные операторы  ,где .

невязка функции - разность между левой и правой частями дифференциального уравнения



Тогда уравнение (1.1) с краевыми условиями (1.2) примет вид:  (4.1)

Решение данного уравнения ищется в виде:

 (4.2)

где  (i=1,2,...,n) - конечная система базисных функций.

Функция y(x) должна удовлетворять краевым условиям, т.е.:



Краевые условия должны выполнятся для любого набора констант Ci, поэтому

 

И исходя из условий курсовой работы () получим:

 

В качестве функции U0(x) выбирается функция



при подстановке ее в полученные выше краевые условия определяются значения констант c и d:



В качестве базисных функций Ui выбираем:



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

При получении точного решения дифференциального уравнения функция .Для приближенного решения невязка отлична от нуля. Для решения близкого к точному, следует подобрать так, чтобы функция  принимала наименьшие значения. Невязка должна быть ортогональна базисным функциям. Таким образом, получается система для определения констант Ci:

 (i=1,2,...,n)

или после подстановки R(x):



что в матричном виде записывается:

 (4.3)

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

Определенные интегралы системы уравнений (4.3) вычислялись по методу Симпсона, общая формула которого имеет вид:



где n- количество разбиений отрезка [a,b] на равные части длины h.

Общая погрешность метода составляет  ,где .

После этого систему уравнений (4.3) разрешается относительно Ci по методу Гаусса (точный метод сложности порядка ).

По полученным значениям Ci и функциям  определяют столбец значений функции y(x) согласно формуле (4.2)

4.3. Метод достижения заданной точности

Достижение заданной точности осуществляется за счет увеличения числа базисных функций, т.е.:

Вычисляем y(x) два раза: один раз с числом базисных функций i(согласно выбранным функциям Ui) другой с i+1, получая при этом значения более точные. Сравниваем значения в первом и во втором случаях. Если максимальная погрешность (норма разности между двумя полученными соответствующими значениями) не превышает заданной точности E, то выбранное число функций Ui можно считать достаточным, полученная функция y(x) удовлетворяет заданной точности и является ответом для данной задачи. В противном случае продолжают увеличивать i, пока не будет достигнута заданная точность.

4.4. Описание результатов

Программа выдает следующие результаты:

Решение краевой задачи для дифференциального уравнения

II-го порядка с точностью Е = 0.00001 методом Галёркина.

x = 1.60 y(x) = -2.000000 eps = 0.0e+00

x = 1.65 y(x) = -2.018974 eps = 4.7e-06

x = 1.70 y(x) = -2.048390 eps = 3.2e-08

x = 1.75 y(x) = -2.090074 eps = 1.8e-06

x = 1.80 y(x) = -2.146268 eps = 3.4e-06

x = 1.85 y(x) = -2.219746 eps = 9.6e-07

x = 1.90 y(x) = -2.313939 eps = 1.6e-06

x = 1.95 y(x) = -2.433100 eps = 4.3e-07

x = 2.00 y(x) = -2.582510 eps = 1.3e-06

x = 2.05 y(x) = -2.768744 eps = 3.2e-06

x = 2.10 y(x) = -3.000000 eps = 0.0e+00

Результаты представлены с шагом 0.05 и максимальной точностью 4,7е-06

В полученной таблице приведено решение заданного дифференциального уравнения 2-ого порядка. В столбце х приведено разбиение отрезка [1,6;2,1] с шагом h=0.05. Погрешность вычисленных значений функции определяется максимальным значением столбца Eps: max_eps, которое должно быть меньше заданной точности E= 0.00001

Приложение 1

Листинг программы “Задача Коши”

#include

#include

#include

#defineN 500

#define epsel 0.00001

double p(double x) {return (-(x));}

double q(double x) {return (-3.1*(x));}

double r(double x) {return (3*(x)+2);}

double a=1.6, b=2.1, A=-2, B=-3;

int NV=10;

void Runge_Kutt ();

int n;

double y[N];

double h, x;

void main()

{

int i, t;

static double y1[N], eps[N];

double max_eps;

for (i=0;iy[i]=0;

n=NV/2;

do

{

for (i=0;i<=n;i++) y1[i]=y[i];

n=n*2;

h=(b-a)/n;

if(n+1>N)

{printf (" ERROR!!!:vihod za predely massiva!");

getch ();

}

Runge_Kutt ();

max_eps=0;

for (i=0;i<=n/2;i++)

{

eps[i]=fabs(y[2*i]-y1[i]);

if (eps[i]>max_eps) max_eps=eps[i];

}

}

while(max_eps>epsel);

clrscr ();

printf("Rezultaty predstavleny s shagom 0.05 ");

printf(" ");

for (i=0;i<=n/2;i=i+n/2/NV)

printf (" x=%.2f y(x)=%+f eps=%.2e ",a+h*2*i,y1[i],eps[i]);

printf (" ");

printf ("Resheno s shagom h=%.3f i max tochnostu E=%.3g metodom Runge Kutta ",h,max_eps);

printf("Chislo razbienii n=%.2d",n);

printf (" Press Any Key to continue...");

getch ()

;

}

void Runge_Kutt ()

{ int i;

static double u[N], v[N], z[N], g[N];

double x,c;

double ku1, ku2,ku3, ku4, kz1, kz2, kz3, kz4;

double tv1, tv2, tv3, tv4, tg1, tg2, tg3, tg4;

z[0]=1;

u[0]=0;

v[0]=A;

g[0]=0;

if(n+1>N)

{printf(" Runge Kutte: small array ");

return;

}

for (i=0;i{x=a+h*i;

ku1=z[i];

kz1=-p(x)*z[i]-q(x)*u[i];

tv1=g[i];

tg1=-p(x)*g[i]-q(x)*v[i]+r(x);

ku2=z[i]+h/2*kz1;

kz2=-p(x+h/2)*(z[i]+h/2*kz1)-q(x+h/2)*(u[i]+h/2*ku1);

tv2=g[i]+h/2*tg1;

tg2=-p(x+h/2)*(g[i]+h/2*tg1)-q(x+h/2)*(v[i]+h/2*tv1)+r(x+h/2);

ku3=z[i]+h/2*kz2;

kz3=-p(x+h/2)*(z[i]+h/2*kz2)-q(x+h/2)*(u[i]+h/2*ku2);

tv3=g[i]+h/2*tg2;

tg3=-p(x+h/2)*(g[i]+h/2*tg2)-q(x+h/2)*(v[i]+h/2*tv2)+r(x+h/2);

ku4=z[i]+h*kz3;

kz4=-p(x+h)*(z[i]+h*kz3)-q(x+h)*(u[i]+h*ku3);

tv4=g[i]+h*tg3;

tg4=-p(x+h)*(g[i]+h*tg3)-q(x+h)*(v[i]+h*tv3)+r(x+h);

u[i+1]=u[i]+(ku1+2*(ku2+ku3)+ku4)*h/6;

z[i+1]=z[i]+(kz1+2*(kz2+kz3)+kz4)*h/6;

v[i+1]=v[i]+(tv1+2*(tv2+tv3)+tv4)*h/6;

g[i+1]=g[i]+(tg1+2*(tg2+tg3)+tg4)*h/6;

}

c=(B-v[n])/u[n];

for (i=0;i<=n;i++) y[i]=u[i]*c+v[i];

}

Приложение 2

Листинг программы “Метод конечных разностей”

#include

#include

#include

#include

#define N 500

#define NV 10

#define EPS 0.00001

#define A -2

#define B -3

int n;

double a=1.6,b=2.1;

double h=0.05,x,y[N];

double eps[N];

double y1[N],max_eps;

void Razn();

double fpx(double);

double fqx(double);

double ffx(double);

double fpx(double x)

{

return (-(x));

}

double fqx(double x)

{

return -3.1*(x);

}

double ffx(double x)

{

return (3*(x)+2);

}

/*double l(double x)

{

return (fqx(x)-2/(h*h))/(1/(h*h)+fpx(x)/(h*2));

}

double m(double x)

{

return (1/(h*h)-fpx(x)/(2*h))/(1/(h*h)+fpx(x)/(h*2));

}

double f(double x)

{

return (ffx(x)/(1/(h*h)+fpx(x)/(h*2));

}*/

main()

{

clrscr();

int i,t;

for(i=0;iy[i]=0;

n=NV/2;

do

{

for(i=0;i<=n;i++) y1[i]=y[i];

n=n*2;

h=(b-a)/n;

if(n+1>N)

{

printf(" Critical Error ");

getch();

return 0;

}

Razn();

max_eps=0;

for(i=0;i<=n/2;i++)

{

eps[i]=fabs(y[2*i]-y1[i]);

if(eps[i]>max_eps) max_eps=eps[i];

}

printf("max eps=%f n=%d ",max_eps,n);

}

while(max_eps>EPS);

clrscr();

printf("Rezultaty podstavlenniye s shagom 0.05 ");

printf(" ");

n=n/2;

h=(b-a)/n;

for(i=0;i<=n;i+=n/NV)

printf("x=%.2f y(x)=%+f eps=%.2e ",a+i*h,y1[i],eps[i]);

printf(" n=%d ",n);

printf("Resheno s shagom h=%f i max. to4nostyu E=%g ",h,max_eps);

printf(" Press any key to continue... ");

getch();

return 0;

}

void Razn()

{

int i;

static double c[N],d[N],u[N],f[N],m[N],th;

y[0]=A;

y[n]=B;

//c[0]=0;

h=(b-a)/n;

for(i=0;i<=n;i++)

{x=a+h*i;

th=1/(h*h)+fpx(x)/(2*h);

u[i]=(fqx(x)-2/(h*h))/th;

m[i]=(1/(h*h)-fpx(x)/(2*h))/th;

f[i]=ffx(x)/th;

}

c[1]=1/u[1];

d[1]=f[1]-m[1]*A;

for(i=1;i{

d[i+1]=f[i+1]-m[i+1]*c[i]*d[i];

c[i+1]=1/(u[i+1]-m[i+1]*c[i]);

}

for(i=n-1;i>=1;i--) {y[i]=c[i]*(d[i]-y[i+1]);}

}

Приложение 3

Листинг программы “Метод Галеркина”

#include

#include

#include

#define P(x) (-(x))

#define Qn(x) (-3.1*(x))

#define R(x) (3*(x)+2)

#define M 11

main()

{

double a=1.6,b=2.1,A=-2,B=-3,eps=0.00001;

double x,y[M][2],u,l,c[M],h,s1,s2,w,norm,n1[M],I[M][M],J[M];

int N=10,Q=10,i,p=0,j,v,s,n=2;

clrscr();

do

{

h=(b-a)/N;n1[0]=0;n1[10]=0;s1=0;s2=0;y[0][p]=A;norm=0;

for(v=1;v<=n;v++)

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

{

for(j=1;j {

x=j*h+a;

u=pow(x-a,i)*(x-b);

l=v*pow(x-a,v-2)*((v-1)*(x-b)+2*(x-a))+P(x)*(v*pow(x-a,v-1)*(x-b)+pow(x-a,v))+Qn(x)*pow(x-a,v)*(x-b);

if(float(j)/2!=j/2) s1+=u*l;

else s2+=u*l;

}

I[i][v]=h/3*(4*s1+2*s2);

s1=0;s2=0;

}

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

{

for(j=1;j {

x=j*h+a;

u=pow(x-a,i)*(x-b);

l=P(x)*(B-A)/(b-a)+Qn(x)*((A*b-B*a)/(b-a)+(B-A)/(b-a)*x);

if(float(j)/2!=j/2) s1+=u*(R(x)-l);

else s2+=u*(R(x)-l);

}

J[i]=h/3*(4*s1+2*s2);

s1=0;s2=0;

}

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

{

w=I[i][s]/I[s][s];

J[i]-=J[s]*w;

for(j=s+1;j<=n;j++) I[i][j]-=I[s][j]*w;

}

c[n]=J[n]/I[n][n];

for(i=n-1;i>0;i--)

{

c[i]=J[i];

for(j=i+1;j<=n;j++) c[i]-=I[i][j]*c[j];

c[i]/=I[i][i];

}

for(j=1;j<=Q;j++)

{

x=j*(b-a)/Q+a;

y[j][p]=(A*b-B*a)/(b-a)+(B-A)/(b-a)*x;

for(i=1;i<=n;i++) y[j][p]=y[j][p]+c[i]*pow(x-a,i)*(x-b);

}

if(p==0) {p=1;n++;norm=eps+1;continue;}

else

{

for(i=1;i {

n1[i]=fabs(y[i][0]-y[i][1]);

if(n1[i]>norm) norm=n1[i];

}

for(i=0;i<=Q;i++) y[i][0]=y[i][1];

}

n++;

}

while (norm>eps);

printf (" ђҐиҐ­ЁҐ Єа Ґў®© § ¤ зЁ ¤«п ¤ЁддҐаҐ­жЁ «м­®Ј® га ў­Ґ­Ёп "

" II-Ј® Ї®ап¤Є  б в®з­®бвмо … = %g ¬Ґв®¤®¬ ѓ «саЄЁ­ . ",eps);

printf (" ?¬Ґп %d га ў­Ґ­Ёп ",n-1);

for(i=0;i<=Q;i++)

printf (" x = %.2f y(x) = %+f eps = %.1e ",i*(b-a)/Q+a,y[i][0],n1[i]);

printf (" ђҐ§г«мв вл ЇаҐ¤бв ў«Ґ­л б и Ј®¬ 0.05 ");

getch();

return 0;

}

Выводы.

При решении ОДУ более рационально использовать метод сведения краевой задачи к задаче Коши или метод конечных разностей.

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

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

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

Литература:

1. Б.П. Демидович и И.А. Марон «Основы вычислительной математики», М. 1970г.