Как известно, метод Ньютона сводится к итерационным вычислениям по следующей формуле (файл р9):
xi+1 = xi +f(xi)/f'(xi).
Реализующая его процедура выглядит довольно просто:
> restart:
NI:=proc(expr, х) local iter;
iter:=x-expr/diff(expr,х);
unapply(iter,x)
end;
NI :=
proc(expr,
x)
local iter, iter:= x - expr/diff(expr,
x); unapply(iter,
x)
end procДля получения итерационной формулы в аналитическом виде здесь используется функция unapply. Теперь, если задать решаемое уравнение, то можно получить искомое аналитическое выражение:
> expr:=sin(х)^2-0.5;
expr := sin(x)² - .5
> F:=NI(expr,x);
Далее, задав начальное приближение для х в виде х=х0, можно получить результаты вычислений для ряда итераций:
> х0:=0.2;
х0:= .2
> to 8 do х0:=F(х0);od;
х0:= 1.382611210 х0:= .117460944 х0:= 2.206529505 х0:= 2.360830634 х0:= 2.356194357 х0:= 2.356194490 х0:= 2.356194490 х0:= 2.356194490
Нетрудно заметить, что, испытав скачок в начале решения, значения х довольно быстро сходятся к конечному результату, дающему корень заданной функции. Последние три итерации дают одно и то же значение х. Заметим, что этот метод дает только одно решение, даже если корней несколько. Вычислить другие корни в таком случае можно, изменив начальное условие.
Можно попробовать с помощью полученной процедуры получить решение и для другой функции:
> expr:=ln(х^2)-0.5;
expr: = ln(x² ) - .5
> F:=NI(expr, х);
F:=x→x-½(ln(x²) - .5)х
>> х0:=0.2;
х0:= 2
> to 8 do x0:=F(x0);od;
x0:= .55718875825 x0:= 1.034437603 x0:= 1.258023119 x0:= 1.283760340 x0:= 1.284025389 x0:= 1.284025417 x0:= 1.284025416 x0:= 1.284025417
Здесь итерационная формула имеет (и вполне естественно) уже другой вид, но сходимость к корню также обеспечивается за несколько итерации.
Возможна и иная форма задания итерационной процедуры с применением оператора дифференцирования D и заданием исходной функции также в виде процедуры (файл р9):
> MI:=proc(f::procedure)
> (х->х)-eval(f)/D(eval(f));
> end;
MI :=
proc(f::procedure) (x→x) -eval(f)/D(eval(f))
end proc> g:=x->x-cos(x);
g := x→x - cos(x)
> SI:=MI(g);
> x0:=0.1;
x0:=.1
> to 6 do x0:=SI(x0) od;
x0:= .9137633858 x0:= .7446642419 x0:= .7390919660 x0:= .7390851333 x0:= .7390851332 x0:= .7390851332
Вообще говоря, в программных процедурах можно использовать любые операторы и функции, присущие Maple-языку, в том числе и те, которые реализуют символьные вычисления. Это открывает широкий простор для разработки новых процедур и функций, обеспечивающих выполнение символьных операций.
10.6.2. Вычисление интеграла по известной формуле
Рассмотрим следующий пример (файл р9):
> Int(e^x*x^n,x)=int(e^x*x^n,x);
∫exxndх = -(-1)(-n)ln(e)(-1-n)(хn(-1)nln(е)nnГ(n)(-х ln(е))(-n) - хn(-1)n ln(е)n е(хln(e)) - хn(-1)nln(е)nn(-х ln(е))(-n)Г(n, -x ln(е)))
Ранние версии системы Maple не брали этот интеграл, поскольку он не имеет аналитического представления через обычные функции. Maple блестяще вычисляет этот «крепкий орешек», но полученное выражение довольно сложно.
Из математики известно, что такой интеграл может быть представлен в следующем виде:
Используя эту формулу, мы можем создать простую процедуру для численного и аналитического вычисления данного интеграла.
> IntExpMonomial:=proc(n::anything,х::name)
local i;
n!*exp(x)*sum(((-1)^(n-i)*х^i)/i!,i=0..n);
end;
IntExpMonomial := proc(n::anything, x::name)
local i;
n! × exp(x) × sum((-1)^(n-i) × x^i/i!,i=-..n)
end proc
Проверим ее в работе:
> IntExpMonomial(3,х);
> IntExpMonomial(5,x);
> IntExpMonomial(n, x);
Результат в аналитическом виде довольно прост для данного интеграла с конкретным значением n. Более того, мы получили несколько иной результат и для n в общем случае. Но точен ли он? Для ответа на этот вопрос продифференцируем полученное выражение:
> diff(%,х);
Результат дифференцирования выглядит куда сложнее, чем вычисленный интеграл. Однако с помощью функции simplify в Maple 9 он упрощается к подынтегральной функции:
> simplify(%);
еx хn
Maple 9.5 выдал более замысловатое выражение:
(-1)n еx (-x)n
Это говорит о том, что задача вычисления заданного интеграла в аналитической форме действительно решена. А что касается громоздкости результатов, так ведь системы, подобные Maple, для того и созданы, чтобы облегчить нам работу с громоздкими вычислениями — в том числе аналитическими.
10.6.3. Вложенные процедуры и интегрирование по частям
Теперь мы подошли к важному моменту, о котором читатель наверняка уже давно догадался — в составляемых пользователям процедурах можно использовать ранее составленные им (или кем-то еще) другие процедуры! Таким образом, Maple-язык позволяет реализовать процедуры, вложенные друг в друга.
Для иллюстрации применения вложенных процедур рассмотрим операцию интегрирования по частям. Пусть нам надо вычислить интеграл
∫ p(x)exdx,
где р(х) — выражение, представляющее полином.
Вначале подготовим процедуру IntExpMonomialR, реализующую вычисление уже рассмотренного ранее интеграла, но рекурсивным способом (файл р9):
> IntExpMonomialR:=proc(n::nonnegint,х::name)
local i;
if n=0 then RETURN(exp(x)) fi;
х^n*ехр(x)-n*IntExpMonomialR(n-1, x);
end;
IntExpMonomialR: = proc(n::nonneg int, x::name)
local i;
if n = 0 then RETURN(exp(x)) end if;
x^n × exp(x) -n × IntExpMonomialR(n - 1, x)
end proc
Проверим ее в работе:
> IntExpMonomialR(4, х);
x4еx - 4x3еx + 12x2ex - 24хеx + 24 еx
> collect(%,exp(х));
(х4 - 4х3 + 12x2 - 24x + 24) еx
Теперь составим процедуру для вычисления по частям нашего интеграла:
> IntExpPolynomial:=proc(р::polynom,x::name) local i, result;
### WARNING: degree(0,x) now returns -infinity
result:=add(coeff(p,x,i)*IntExpMonomialR(i,x),i=0..degree(p, x));
collect(result,exp(x));
end;
IntExpPolynomial: = proc(p::polynom, x::name)
local i, result;
result:= add(coeff(p, x, i) × IntExpMonomialR(i, x), i = 0..degree(p, x)); collect(result, exp(x))
end proc
В этой процедуре имеется обращение к ранее составленной процедуре IntExpMonomialR. Обратите внимание на то, что в процедуре введено предупреждение об определенных проблемах, связанных с использованием функции degree (сообщение начинается с символов ###). Тем не менее, процедура работает, в чем убеждают, по крайней мере, следующие примеры:
> р:=(х^2+1)*(1-3*х);
р := (х² + 1)(1 - 3х)
> expand(р);
х² - 3х³ + 1 - 3x
> int(р*ехр(х),х);
-ех(-24 + 23х - 10х² + 3х³)
> IntExpPolynomial(р,х);
(24 - 23х + 10х² - 3х³)еx
> IntExpPolynomial(р,х);
(24 - 23х + 10х² - 3х³)еx
10.7. Дополнительные возможности Maple-языка
10.7.1. Переназначение определений
В большинстве случаев Maple-язык использует достаточно длинные идентификаторы для своих определений, например функций. Однако с помощью функции alias можно изменить любое определение на другое, если оно кажется пользователю более удобным. Функция alias записывается в виде
alias(e1, е2, ..., eN)
где е1, е2, …, eN — ноль или более равенств.
Эта функция возвращает список переназначений и осуществляет сами переназначения. Например, для замены имени функции BesselJ на более короткое имя BJ достаточно параметром функции alias записать BJ=BesselJ:
> alias(ВJ=BesselJ);
BJ, Fx
> [BJ(0,1.),BesselJ(0,1.)];
[.7651976866, .7651976866]
Можно также переназначить функцию пользователя:
> alias(Fx=F(x));
