Некоторые замечательные функции
6.1. Функция Бесселя
Для
Легко видеть, что
Для определения слагаемых применим рекуррентную формулу
В качестве начальных значений суммы и слагаемого использовать член ряда при k = 0, равный
Для вычисления значения функции при k = 0, необходимо использовать функции вычисления степени c основанием
Функция вычисления степени.
Function
Extent(x : real; n : integer) : real;
var
i : integer;
E : real;
begin
E := 1;
if n = 0 then E := 1
else for
i := 1 to n do E := E*x/2;
Extent := E
end;
Функция вычисления факториала.
Function Fakt(n : integer) : longint;
var
i : integer;
F : longint;
begin
F := 1;
if n = 0 then F := 1
else for
i := 1 to n do F := F*i;
Fakt := F
end;
А теперь составим функцию вычисления всей суммы
Function J(x, eps : real; n : integer) : real;
var
y, jj : real;
k : integer;
begin
k := 0; y := Extent(x, n)/Fakt(n); jj := 0;
repeat
jj := jj + y;
k := k + 1;
y := -y*x*x/(4*k*(n + k))
until abs(y) < eps;
J := jj
end;
Программа
Program Bessel;
uses WinCrt;
var
n : integer;
x, eps : real;
{----------------------------------------------------------------------------------------}
Function
t(eps : real) : integer;
var
k : integer;
begin
k := -1;
repeat
eps := eps*10;
k := k + 1
until eps > 1;
t := k
end;
{----------------------------------------------------------------------------------------}
Function Extent(x : real; n : integer) : real;
var
i : integer;
E : real;
begin
E := 1;
if n = 0 then E := 1
else for
i := 1 to n do E := E*x/2;
Extent := E
end;
{----------------------------------------------------------------------------------------}
Function Fakt(n : integer) : longint;
var
i : integer;
F : longint;
begin
F := 1;
if n = 0 then F := 1
else for
i := 1 to n do F := F*i;
Fakt := F
end;
{----------------------------------------------------------------------------------------}
Function J(x, eps : real; n : integer) : real;
var
y, jj : real;
k : integer;
begin
k := 0; y := Extent(x, n)/Fakt(n); jj := 0;
repeat
jj := jj + y;
k := k + 1;
y := -y*x*x/(4*k*(n + k))
until abs(y) < eps;
J := jj
end;
{----------------------------------------------------------------------------------------}
begin
write('Введите значение x '); readln(x);
write('Введите значение n '); readln(n);
write('Введите точность вычисления '); readln(eps);
writeln('Значение функции Бесселя равно ', J(x, eps, n):6:t(eps))
end.
6.2. Гамма-функция Эйлера
Рассмотрим (вместе с Эйлером) бесконечное произведение
считая x отличным от нуля и от всех целых отрицательных чисел.
Легко представить его общий множитель так:
отсюда вытекает, что наше произведение (абсолютно) сходится.
Определяемая им функция
Используя формулу
Теперь надо оставить функцию для вычисления степени действительного положительного аргумента a с любым действительным показателем x. Для этого достаточно воспользоваться соотношением:
Function
Extent_real(a, x : real) : real;
begin
Extent_real := exp(x*ln(a))
end;
Теперь можно составить основную функцию для вычисления гамма-функции.
Function G(x, eps : real) : real;
var
n : longint;
g1, gg : real;
begin
n := 1; gg := Extent_real((n + 1)/n, x)/(x*(x + n));
repeat
n := n + 1;
gg := gg*n*Extent_real((n + 1)/n, x)/(x + n);
n := n + 1;
g1 := gg*n*Extent_real((n + 1)/n, x)/(x + n)
until abs(g1 - gg) < eps;
G := gg
end;
Программа
Program
Gamma_function;
uses WinCrt;
var
x, eps : real;
{----------------------------------------------------------------------------------------}
{ Функция вычисления порядка - количества знаков после запятой }
Function t(eps : real) : integer;
var
k : integer;
begin
k := -1;
repeat
eps := eps*10;
k := k + 1
until eps > 1;
t := k
end;
{----------------------------------------------------------------------------------------}
Function Extent_real(a, x : real) : real;
begin
Extent_real := exp(x*ln(a))
end;
{----------------------------------------------------------------------------------------}
Function G(x, eps : real) : real;
var
n : longint;
g1, gg : real;
begin
n := 1; gg := Extent_real((n + 1)/n, x)/(x*(x + n));
repeat
n := n + 1;
gg := gg*n*Extent_real((n + 1)/n, x)/(x + n);
n := n + 1;
g1 := gg*n*Extent_real((n + 1)/n, x)/(x + n)
until abs(g1 - gg) < eps;
G := gg
end;
{----------------------------------------------------------------------------------------}
begin
write('Введите значение аргумента '); readln(x);
write('Введите точность вычисления '); readln(eps);
writeln('Значение гамма-функции равно ', G(x, eps):6:t(eps))
end.
Так как n- е частичное произведение имеет вид
Написав аналогичную формулу для
и мы приходим к простому и важному соотношению:
Если вы выполняли предыдущую программу для различных значений аргумента, то смогли убедиться, что точность вычисления невелика, даже при достаточно большом значении n.
Формула
Для этого составим следующую функцию:
Function
Gamma(x, eps : real) : real;
begin
x := x - 1;
Gamma := x*G(x, eps)
end;
Program
Gamma_function;
uses WinCrt;
var
x, eps, gg : real;
{----------------------------------------------------------------------------------------}
{ Функция вычисления порядка - количества знаков после запятой }
Function
t(eps : real) : integer;
var
k : integer;
begin
k := -1;
repeat
eps := eps*10;
k := k + 1
until eps > 1;
t := k
end;
{----------------------------------------------------------------------------------------}
Function Extent_real(a, x : real) : real;
begin
Extent_real := exp(x*ln(a))
end;
{----------------------------------------------------------------------------------------}
Function G(x, eps : real) : real;
var
n : longint;
g1, gg : real;
begin
n := 1; gg := Extent_real((n + 1)/n, x)/(x*(x + n));
repeat
n := n + 1;
gg := gg*n*Extent_real((n + 1)/n, x)/(x + n);
n := n + 1;
g1 := gg*n*Extent_real((n + 1)/n, x)/(x + n)
until abs(g1 - gg) < eps;
G := gg
end;
{----------------------------------------------------------------------------------------}
Function Gamma(x, eps : real) : real;
begin
x := x - 1;
Gamma := x*G(x, eps)
end;
{----------------------------------------------------------------------------------------}
begin
write('Введите значение аргумента '); readln(x);
write('Введите точность вычисления '); readln(eps);
if (x < 2) and (x > 1) then gg := Gamma(x, eps)
else gg := G(x, eps);
writeln('Значение гамма-функции равно ', gg:6:t(eps))
end.
Если положить x равным натуральному числу m, то получим рекуррентную формулу
Так как
Эта формула дает еще одну возможность уточнить значения функции для натуральных значений аргумента. Для вычисления по этой формуле составим еще одну функцию.
Function G_natural(m : integer) : longint;
var
i : integer;
g : longint;
begin
g := 1;
for i := 2 to m do
g := g*i;
G_natural := g
end;
Программа
изменится и станет такой:
Program
Gamma_function;
uses WinCrt;
var
x, eps, gg : real;
{----------------------------------------------------------------------------------------}
{ Функция вычисления порядка - количества знаков после запятой }
Function
t(eps : real) : integer;
var
k : integer;
begin
k := -1;
repeat
eps := eps*10;
k := k + 1
until eps > 1;
t := k
end;
{----------------------------------------------------------------------------------------}
Function Extent_real(a, x : real) : real;
begin
Extent_real := exp(x*ln(a))
end;
{----------------------------------------------------------------------------------------}
Function G(x, eps : real) : real;
var
n : longint;
g1, gg : real;
begin
n := 1; gg := Extent_real((n + 1)/n, x)/(x*(x + n));
repeat
n := n + 1;
gg := gg*n*Extent_real((n + 1)/n, x)/(x + n);
n := n + 1;
g1 := gg*n*Extent_real((n + 1)/n, x)/(x + n)
until abs(g1 - gg) < eps;
G := gg
end;
{----------------------------------------------------------------------------------------}
Function Gamma(x, eps : real) : real;
begin
x := x - 1;
Gamma := x*G(x, eps)
end;
{----------------------------------------------------------------------------------------}
Function G_natural(m : integer) : longint;
var
i : integer;
g : longint;
begin
g := 1;
for i := 2 to m do
g := g*i;
G_natural := g
end;
{----------------------------------------------------------------------------------------}
begin
write('Введите значение аргумента '); readln(x);
write('Введите точность вычисления '); readln(eps);
if x = trunc(x)
then gg := G_natural(trunc(x))
else
if (x < 2) and (x > 1)
then gg := Gamma(x, eps)
else gg := G(x, eps);
writeln('Значение гамма-функции равно ', gg:6:t(eps));
end.
Еще одну важную формулу для функции
Мы находим:
Это - формула Вейерштрасса.
Библиотека часто встречающихся процедур и функций
55. Вычисление числа e.
Procedure number_e(eps : real; var
e : real);
var
n : integer;
u : real;
begin
e := 0; u := 1; n := 1;
repeat
e := e + u;
u := u/n; n := n + 1
until 3*u <= eps;
end;
56. Вычисление корней любой степени из произвольного числа.
Procedure Radical(n : integer ; x, eps : real; var b : real);
var
n : integer;
z, m, u : real;
begin
b := 1; u := 1; n := 1;
repeat
u := (m - n + 1)*x*u/n; b := b + u; n := n+1
until abs(u) <= eps;
end;
57. Процедура вычисления бесконечного произведения
Procedure
Multiplication(eps : real; var Mult : real);
var
n : longint;
Mult1 : real;
begin
n := 2; Mult1 := 1;
repeat
Mult1 := Mult1*(1 - 1/sqr(n));
n := n + 1;
Mult := Mult1*(1 - 1/sqr(n))
until abs(Mult - Mult1) < eps
end;
58. Процедура вычисления числа
Procedure Wallis(eps : real; var Mult : real);
var
n : longint;
Mult1 : real;
begin
n := 1; Mult := 1;
repeat
Mult := Mult*(4*sqr(n)/(4*sqr(n)-1));
n := n + 1;
Mult1 := 4*sqr(n)/(4*sqr(n)-1)
until Mult1 < eps
end;
59. Процедур вычисления эллиптического интеграла 1-го рода через бесконечное произведение.
Procedure Elliptic(k, eps : real; var Kk : real);
var
Kk1 : real;
begin
Kk1 := k;
repeat
k := (1 - sqrt(1 - sqr(k)))/(1 + sqrt(1 - sqr(k)));
Kk1 := Kk1*(1 + k);
k := (1 - sqrt(1 - sqr(k)))/(1 + sqrt(1 - sqr(k)));
Kk := Kk1*(1 + k);
until abs(Kk1 - Kk) < eps;
Kk := Kk*Pi/2
end;
60. Рекуррентная функция вычисления интеграла вероятностей.
Function
FF(x : real) : real;
var
n : integer;
u, I : real;
begin
if x >= 5
then FF := 1
else if x <= -5
then FF := -1
else
begin
u := x; n := 0; I := 0;
repeat
I := I + u;
n := n + 1;
u := -u*(x*x*(2*n - 1)/(2*n*(2*n + 1)))
until abs(u) < 0.00001;
FF := 2*I/sqrt(2*Pi)
end
end;
61. Процедура вычисления интеграла
Procedure
Integral(x, eps : real; var I : real);
var
n : integer;
u : real;
begin
u := x; n := 1; I := 0;
repeat
I := I + u;
n := n + 1;
u := -(u*x*x*(2*n - 3))/((2*n - 2)*sqr(2*n - 1))
until abs(u) < eps
end;
62. Процедура вычисления эллиптического интеграла 2-го рода с помощью интегрирования рядов.
Procedure Elliptic2(k, eps : real; var Ek : real);
var
n : integer;
u : real;
begin
u := k*k/4; n := 1; Ek := 0;
repeat
Ek := Ek + u;
n := n + 1;
u := (u*k*k*(2*n - 1)*(2*n - 3))/(4*n*n);
until abs(u) < eps;
Ek := Pi*(1 - Ek)/2
end;
63. Функция Бесселя:
Function
J(x, eps : real; n : integer) : real;
var
y, jj : real;
k : integer;
begin
k := 0; y := Extent(x, n)/Fakt(n); jj := 0;
repeat
jj := jj + y;
k := k + 1;
y := -y*x*x/(4*k*(n + k))
until abs(y) < eps;
J := jj
end;
64. Функция вычисления степени положительного действительного числа с произвольным действительным показателем.
Function
Extent_real(a, x : real) : real;
begin
Extent_real := exp(x*ln(a))
end;
65. Функции для вычисления гамма-функции.
Для произвольного действительного аргумента x, отличного от нуля и от всех целых отрицательных чисел.
Function G(x, eps : real) : real;
var
n : longint;
g1, gg : real;
begin
n := 1; gg := Extent_real((n + 1)/n, x)/(x*(x + n));
repeat
n := n + 1;
gg := gg*n*Extent_real((n + 1)/n, x)/(x + n);
n := n + 1;
g1 := gg*n*Extent_real((n + 1)/n, x)/(x + n)
until abs(g1 - gg) < eps;
G := gg
end;
Вычисление гамма-функции для значений аргумента из промежутка (1, 2)
Function Gamma(x, eps : real) : real;
begin
x := x - 1;
Gamma := x*G(x, eps)
end;
Вычисление гамма-функции натурального аргумента.
Function
G_natural(m : integer) : longint;
var
i : integer;
g : longint;
begin
g := 1;
for i := 2 to m do
g := g*i;
G_natural := g
end;
Упражнения
194. Вычислить сумму членов рядов, заданных формулами общего члена. Число членов ряда задается пользователем.
195. Составьте функции для нахождения с указанной точностью числа
a)
б) ряд Леонардо Эйлера:
196. Составьте функции для вычисления с указанной точностью значений тригонометрических функций путем вычисления следующих рядов:
а)
Для вычисления значения члена ряда использовать рекуррентную формулу
б)
Для определения значения члена ряда использовать формулу
197. Дано действительное x. Вычислить приближенное значение бесконечной суммы:
Нужное приближение считается полученным, если вычислена сумма нескольких первых слагаемых, и очередное слагаемое оказалось по модулю меньше 0.00001.
198. Вычислить приближенное значение бесконечной суммы (справа от каждой суммы дается ее точное значение, с которым можно сравнить полученный результат):
а)
б)
в)
г)
Нужное приближение считается полученным, если вычислена сумма нескольких первых слагаемых, и очередное слагаемое оказалось по модулю меньше данного положительного числа eps.
199. Дано: натуральное n, действительное x. Вычислить:
а)
в)
200. Дано: действительные a, h, натуральное n. Вычислить
201. Дано: натуральное n, действительное x. Вычислить:
202. Вычислить:
а)
в)
203. Вычислить сумму членов ряда
с точностью до члена ряда меньшего 0.0000001.
204. Установить, сходятся ли следующие ряды:
а)
205. Составить программу вычисления бесконечного произведения:
Частичное произведение имеет вид
где C - эйлерова постоянная, а
206. Эйлером были найдены следующие разложения тригонометрических функций в бесконечные произведения:
1)
2)
Составьте программы вычисления этих произведений и вычислите значения sinx и cosx с заданной точностью.
207. Дано натуральное n. Получить:
а)
в)
г)
д)
208. Дано натуральное n. Получить
209. Вычислить интеграл
Ответы
К заданию 1
Program
Task1; {Вычисление sinx с помощью ряда}
uses WinCrt;
var
n, k : integer;
x, e, eps, sin, u : real;
{----------------------------------------------------------------------------------------}
{ Функция вычисления порядка - кол-во знаков после запятой }
Function t(eps : real) : integer;
var
k : integer;
begin
k := -1;
repeat
eps := eps*10;
k := k + 1
until eps > 1;
t := k
end;
{----------------------------------------------------------------------------------------}
begin
write('Задайте точность вычисления '); readln(eps);
write('Введите значение аргумента в радианах '); readln(x);
u := x; n := 2;
sin := x;
repeat
u := -u*x*x/((2*n-1)*(2*n-2));
sin := sin + u;
n := n + 1
until abs(u) < eps;
write('Значение sin( ', x:1:t(eps), ' ) равно ', sin:3:t(eps));
writeln(' с точностью до ', eps:1:t(eps))
end.
К
заданию 2
Program Task2;
uses WinCrt;
var
n : integer;
p, u, eps : real;
{----------------------------------------------------------------------------------------}
{ Функция вычисления порядка - кол-во знаков после запятой }
Function
t(eps : real) : integer;
var
k : integer;
begin
k := -1;
repeat
eps := eps*10;
k := k + 1
until eps > 1;
t := k
end;
{----------------------------------------------------------------------------------------}
begin
write(' Укажите точность вычисления числа пи '); readln(eps);
u := 1; p := 0; n := 1;
repeat
p := p + u;
n := n + 1;
u := -u*(2*n - 3)/(2*n - 1)
until abs(u) < eps;
write('Число Пи равно ', 4*p:1:t(eps), ' с точностью до ');
writeln(eps:1:t(eps))
end.
К
заданию 3
Program Task3;
uses WinCrt;
var
k : integer;
a, u, z, z1, e, eps : real;
{----------------------------------------------------------------------------------------}
{ Функция вычисления порядка - кол-во знаков после запятой }
Function t(eps : real) : integer;
var
k : integer;
begin
k := -1;
repeat
eps := eps*10;
k := k + 1
until eps > 1;
t := k
end;
{----------------------------------------------------------------------------------------}
begin
write('Введите значение a>1 '); readln(a);
write('Задайте точность eps '); readln(eps);
u := 1; z := 0;
repeat
u := u*a;
z := z + 1/(1 + u);
z1 := z + 1/(1 + u*a)
until abs(z1 - z) < eps;
write('Сумма ряда равна ', z1:1:t(eps));
writeln(' с точностью до ', eps:1:t(eps))
end.
К
заданию 4
Program Task4;
uses WinCrt;
n : integer;
var
a : longint;
{----------------------------------------------------------------------------------------}
Function c(a : integer) : longint;
var
s, z : integer;
begin
z := 0;
repeat
s := a mod 10;
z := z + s;
a := a div 10
until s = 0;
c := z
end;
{----------------------------------------------------------------------------------------}
begin
a := 7;
n := 1;
while n <= 1000 do
begin
a := c(a*a);
n := n + 1
end;
writeln(n - 1, '- й член последовательности равен ', a)
end.
