Задачи высшей математики с Maxima

Если $$z = x + iy = r  (\cos (j) + i \sin (j) )$$, то значения функции $$f(z) = exp(z)$$ вычисляются по формуле $$f(z) = e^z = e^{x+iy} =  e^x (\cos (y) + i \sin (y))$$.

(%i1)	eq:x^2-4*x+13=0;

(%i2)	solve(eq,x);

(%i3)	x1:%o2[1]$ x2:%o2[2];

(%i4)	print(x1,x2);

(%i1)	solve(x^3=1,x);

(%i2)	solve(x^5=1,x);

(%i1)	carg (1);

(%i2)	carg (1 + %i);

(%i3)	carg (exp (%i)),numer;

(%i4)	carg (exp (%pi * %i));

(%i5)	carg (exp (3/2 * %pi * %i));

(%i1)	demoivre:true;

(%i2)	demoivre (exp (3+3/2 * %pi * %i));

(%i3)	demoivre (exp (%pi+3/2 * %pi * %i));

(%i1)	declare ([aa, bb], real, cc, complex, ii, imaginary);

(%i2)	conjugate (aa + bb*%i);

(%i3)	conjugate (ii);

(%i1)	a:1+%i;

(%i2)	plog(a);

(%i1)	a:1+%i;

(%i2)	polarform(a);

(%i3)	rectform(%);

(%i1)	a:matrix([1,2,3],[4,5,6]);

(%i1)	a:matrix([1,2,3]); transpose(a);

(%i2)	c:b*2;

(%i4)	 b+c;

(%i5)	a+b;

fullmap: arguments must have same formal structure.
– an error. To debug this try: debugmode(true);

(%i6)	f:a.b;

(%i1)	a:matrix([1,2],[3,4]);
	b:invert(a);
	b.a;

(%i4)	determinant(a);

(%i6)	charpoly(a,x);

(%i7)	ratsimp(%);

(%i1)	a:matrix([1,1,1],[2,2,2],[3,3,3]);

(%i2)	eigenvalues(a);

(%i3)	eigenvectors(a);

(%i1)	load("eigen");

(%i2)	x:matrix([1,2,3],[4,5,6]);

(%i3)	y:gramschmidt(x);

(%i4)	ratsimp(%[1].%[2]);

(%i1)	a:matrix([1,2,3],[4,5,x],[6,7,y]);

(%i2)	b:echelon(a);

(%i1)	a:matrix([1,2,3,4],[2,5,6,9]);

(%i1)	a:matrix([1,2,3,4],[2,5,6,9]);

(%i2)	rank(a);

(%i3)	b:matrix([1,1],[2,2],[3,3],[4,5]);

(%i4)	rank(b);

(%i1)	A: matrix( [1, 2, 2],[ -1, -1, 3], [2, 5, 0]);

(%i2)	B:matrix([10, 0],[-2, 5], [1, 4]);

(%i3)	determinant(A);

(%i4)	A1:invert(A); x:A1.B;

(%i6)	A.x-B;

(%i10)	C:matrix([10,0,-2],[5,1,4]); x:C.A1; x.A-C;

(%i1)	ex1:x+y+z+t=6;
	ex2:2*x-2*y+z+3*t=2;
	ex3:3*x-y+2*z-t=8;
	linsolve([ex1,ex2,ex3],[x,y,z,t]);

(%i5)	ev(%),%r1=1;

(%i6)	e1: 2*x*(1 - a1) - 2*(x - 1)*a2; e2: a2 - a1;
	e3: a1*(-y - x^2 + 1); e4: a2*(y - (x - 1)^2);

(%i10)	algsys ([e1, e2, e3, e4], [x, y, a1, a2]);

(%i11)	realroots (2 - x + x^5, 5e-6);

(%i12)	float(%);

(%i13)	ev(2-x+x^5,%[1]);

(%i1)	eqn:x^4+1; soln:allroots (eqn);

(%i1)	p: x^10 - 2*x^4 + 1/2$ nroots (p, -6, 9);

(%i1)	eqn:x^2+x+1=(x-1)^3;

(%i2)	lhs(eqn);

(%i3)	rhs(eqn);

(%i1)	expr1: 2*x^2 + y*x + z;
	expr2: 3*x + 5*y - z - 1;
	expr3: z^2 + x - y^2 + 5;

(%i4)	eliminate ([expr3, expr2, expr1], [y, z]);

(%i1)	limit(sin(x)/x, x, 0);

(%i2)	limit (exp(x), x, inf);

(%i3)	limit (exp(x), x, minf);

(%i3)	limit(1/x, x, 0, minus);

(%i4)	limit(1/x, x, 0, plus);

(%i8)	limit(x^(1/3), x, 8);

(%i9)	limit(x^(1/3), x, -8);

(%i10)	y(x):=(x^3-3*x-2)/(x^2-x-2)^2; limit(y(x), x,-1);

(%i16)	limit((x^2-4)/(x^2-3*x+2), x, 2);

(%i17)	factor((x^2-4)/(x^2-3*x+2));

(%i18)	limit(%,x,2);

(%i1)	limit((sqrt(x)-1)/(x-1), x, 1);

(%i2)	factor((sqrt(x)-1)/(x-1));

(%i3)	radcan(%);

(%i4)	limit(%, x, 1);

(%i1)	limit(sin(5*x)/x, x, 0);

(%i4)	limit((1-cos(2*x))/x^2, x, 0);

(%i5)	limit(log(1+x)/x, x, 0);

(%i6)	limit((a^x-1)/x, x, 0);

(%i7)	limit((x/(2+x))^(3*x), x, inf);

(%i26)	f(x):=sin(3*x)*sin(5*x)$ g(x):=(x-x^3)^2$

(%i28)	limit(f(x)/g(x), x, 0);

(%i1)	tlimswitch=true;

(%i1)	f(x):=(tan(x)-sin(x))/(x-sin(x));

(%i2)	tlimit(f(x),x,0);

3.3.3 Непрерывные функции

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

3.3.3.1 Непрерывность функции в точке

Дадим два определения понятия непрерывности функции в точке.

(рис 3.1) Определение непрерывной функции

Определение 1. Функция $$f(x)$$ называется непрерывной в точке $$a$$, если она удовлетворяет трем условиям: 1) $$f(x)$$ определена в некоторой окрестности точки $$x = a$$, 2) существует конечный предел l$$\displaystyle{\lim_{x \to a} f(x)}$$ этот предел равен значению функции $$f(x)$$ в точке $$a$$, т.е. $$\displaystyle{\lim_{x\to a}f(x)}=f(a)$$. Очевидно, что непрерывность функции в данной точке выражается непрерывностью ее графика при прохождении данной точки.

Рассмотрим второе определение непрерывности функции в точке.

Придадим аргументу $$a$$ приращение $$\Delta x\ne 0$$. Тогда функция $$y = f(x)$$ получит приращение $$\Delta y$$, определяемое как разность наращенного и исходного значения функции: $$\Delta y=f(a+\Delta x)-f(a)$$ (см. рис. 3.1).

Определение 2. Функция $$y = f(x)$$ называется непрерывной в точке $$a$$, если она определена в некоторой окрестности точки $$x = a$$, и приращение ее $$\Delta y$$ в этой точке, соответствующее приращению $$\Delta x$$, стремится к нулю при стремлении $$\Delta x$$ к нулю: $$\displaystyle{\lim_{\Delta x\to 0}\Delta y}=0$$.

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

Пример исследования непрерывности функции с Maxima:

Функция $$f\left( x\right) :=\frac{1}{1+exp\left( \frac{1}{1-x}\right) }$$ имеет возможную точку разрыва при $$x = 1$$. Сопоставим пределы данной функции при стремлении $$x$$ к 1 слева и справа:

(%i16)	f(x):=1/(1+exp(1/(1-x)));

$$f\left( x\right) :=\frac{1}{1+exp\left( \frac{1}{1-x}\right) }\leqno{(\%o16) }$$

(%i17)	limit(f(x),x,1,plus);

$$1\leqno{(\%o17) }$$

(%i18)	limit(f(x),x,1,minus);

$$0\leqno{(\%o18) }$$

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

3.3.3.2 Свойства непрерывных функций

Если функции $$f(x)$$ и $$g(x)$$ непрерывны в точке $$a$$, то их сумма $$f(x) + g(x)$$, произведение $$f(x)g(x)$$, и частное $$\frac{f(x)}{g(x)}$$ (при условии, что $$g(a)\ne 0$$) являются функциями, непрерывными в точке $$a$$.

Если функция $$y = f(x)$$ непрерывна в точке $$a$$ и $$f(a)>0$$, то существует такая окрестность точки a, в которой $$f(x)>0$$.

Если функция $$y = f(u)$$ непрерывна в точке $$u_0$$, а функция $$u=\psi(x)$$ непрерывна в точке $$u_0=\psi(x_0)$$, то сложная функция $$y=f[\psi(x)]$$ непрерывна в точке $$x_0$$.

Свойство 3 может быть записано в виде: $$\lim_{x\to x_0}f[\psi(x)]=f\left[\lim_{x\to x_0}\psi(x)\right],$$ т.е. под знаком непрерывной функции можно переходить к пределу.

Функция $$y = f(x)$$ называется непрерывной на промежутке $$X$$, если она непрерывна в каждой точке этого промежутка. Можно доказать, что все элементарные функции непрерывны в области их определения.

3.3.3.3 Точки разрыва функций и их классификация

Точка $$a$$, принадлежащая области определения функции или являющаяся граничной для этой области, называется точкой разрыва функции $$f(x)$$, если в этой точке нарушается условие непрерывности функции.

Если существуют конечные пределы $$f(a-0)=\displaystyle{\lim_{x\to a-0}f(x)}\; \text{и}\; f(a+0)=\displaystyle{\lim_{x\to a+0}f(x)},$$ причем не все три числа $$f(a), f(a - 0), f(a + 0)$$ равны между собой, то точка $$a$$ называется точкой разрыва 1 рода (существуют конечные односторонние пределы функции слева и справа, не равные друг другу).

Точки разрыва 1 рода подразделяются, в свою очередь, на точки устранимого разрыва (когда $$f(a - 0) = f(a + 0) \ne f(a)$$, т.е. когда левый и правый пределы функции $$f(x)$$ в точке a равны между собой, но не равны значению функции $$f(x)$$ в этой точке) и на точки скачка (когда $$f(a-0) \ne f(a+0)$$, т.е. когда левый и правый пределы функции в точке $$a$$ различны); в последнем случае разность $$f(a + 0) - f(a - 0)$$ называется скачком функции $$f(x)$$ в точке $$a$$.

Точки разрыва, не являющиеся точками разрыва 1 рода, называются точками разрыва 2 рода. В точках разрыва 2 рода не существует хотя бы один из односторонних пределов.

Рассмотрим предыдущий пример. Функция $$f\left( x\right) :=\frac{1}{1+exp\left( \frac{1}{1-x}\right) }$$ имеет точку разрыва при $$x = 1$$.

(рис 3.2) Разрыв исследуемой функции

Так как пределы $$\lim_{x\to {0-0}} f \left( x\right)$$ и $$\lim_{x\to {0+0}} f \left( x\right)$$ не совпадают, но оба конечны, делаем вывод о наличии точки разрыва первого рода при $$x = 1$$.

Графическую иллюстрацию получаем при помощи wxMaxima (см. рис. 3.2).

3.3.4 Дифференцирование с помощью пакета Maxima

Пакет Maxima предоставляет мощные средства для дифференцирования функций и вычисления дифференциалов. Для вычисления простейшей производной следует в командном окне после приглашения Maxima ввести команду следующего вида: diff(<функция>, <переменная>); где <функция> –– выражение, задающее функцию (не обязательно одной переменной), например x^2+2*x+1; <переменная> –– имя переменной, по которой будет вестись дифференцирование, например $$x$$.

Примером вычисления производной может служить такая команда: diff(x^2+2*x+1,x);.

С помощью команды $$diff$$ можно вычислять производные высших порядков. При этом команда имеет следующий формат: $$diff$$(<функция>,<переменная>,<порядок>); где <порядок> — порядок вычисляемой производной.

В решениях некоторых примеров этой главы с помощью Maxima будут использованы дополнительные команды Maxima:

$$ratsim$$(<выражение >), $$radcan$$(<выражение >) — упрощение алгебраического выражения.

$$trigsim$$(<выражение >), $$trigexpand$$(<выражение >) — упрощение или подстановка тригонометрического выражения.

$$factor$$(<выражение>); –– разложить <выражение> на множители.

$$at$$(<выражение>,<old>=<new>); — подставить выражение <new> на место <old> в <выражении>.

<переменная>:$$solve$$(<выр1>=<значение>,<выр2>); —- присвоить <переменной> значение выражения <выр2>, полученное разрешением уравнения <выр1>(<выр2>)=<значение>.

$$taylor(f(x),x,x0,n);$$ — разложить функцию $$f(x)$$ по формуле Тейлора с центром в точке $$x0$$ до порядка n включительно.

3.3.4.1 Вычисление производных и дифференциалов

Для вычисления производной функции используется функция $$diff$$, для вычисления производных различного порядка удобно создать пользовательскую функцию (в примере ниже — $$f(x)$$):

(%i3)	f(x):=sin(9*x^2);

$$f\left( x\right) :=sin\left( 9\,{x}^{2}\right) \leqno{(\%o3) }$$

(%i4)	d1:diff(f(x),x,1);

$$18\,x\,cos\left( 9\,{x}^{2}\right) \leqno{(\%o4) }$$

(%i5)	d2:diff(f(x),x,2);

$$18\,cos\left( 9\,{x}^{2}\right) -324\,{x}^{2}\,sin\left( 9\,{x}^{2}\right) \leqno{(\%o5) }$$

(%i6)	d3:diff(f(x),x,3);

$$-972\,x\,sin\left( 9\,{x}^{2}\right) -5832\,{x}^{3}\,cos\left( 9\,{x}^{2}\right) \leqno{(\%o6) }$$

Пример вычисления дифференциала ($$del(x)$$ равноценно $$dx$$, не указана явно переменная дифференцирования):

(%i8)	diff(log(x));

$$\frac{del\left( x\right) }{x}\leqno{(\%o8) }$$

Аналогичный подход применим и для функции нескольких переменных. Функция $$diff$$ с единственным аргументом — дифференцируемой функцией — возвращает полный дифференциал.

Пример:

(%i9)	diff(exp(x*y));

$$x\,{e}^{x\,y}\,del\left( y\right) +y\,{e}^{x\,y}\,del\left( x\right) \leqno{(\%o9) }$$

Пример:

(%i10)	diff(exp(x*y*z));

$$x\,y\,{e}^{x\,y\,z}\,del\left( z\right) +x\,z\,{e}^{x\,y\,z}\,del\left( y\right) +y\,z\,{e}^{x\,y\,z}\,del\left( x\right) \leqno{(\%o10) }$$

Если указать апостроф перед символом $$diff$$, то производная не вычисляется и упрощение, обычно предусмотренное по умолчанию, не осуществляется.

Пример:

Создаём функцию $$f(x,z)$$:

(%i18)	f(x,z):=x^2*z+z^2*x;

$$f\left( x,z\right) :={x}^{2}\,z+{z}^{2}\,x\leqno{(\%o18) }$$

Вычисляем дифференциальное выражение:

(%i19)	diff (f(x,z), x, 2) + diff (f(x,z), z, 3) +
	diff (f(x,z), x) * x^2;

$${x}^{2}\,\left( {z}^{2}+2\,x\,z\right) +2\,z\leqno{(\%o19)}$$

Производим формальное дифференцирование, не вычисляя непосредственно результат:

(%i20)	'diff (f(x,z), x, 2) + 'diff (f(x,z), z, 3) +
	'diff (f(x,z), x) * x^2;

$$\frac{{d}^{3}}{d\,{z}^{3}}\,\left( x\,{z}^{2}+{x}^{2}\,z\right) +\frac{{d}^{2}}{d\,{x}^{2}}\,\left( x\,{z}^{2}+{x}^{2}\,z\right) +{x}^{2}\,\left( \frac{d}{d\,x}\,\left( x\,{z}^{2}+{x}^{2}\,z\right) \right) \leqno{(\%o20) }$$

3.4 Экстремумы функций

3.4.1 Отыскание максимумов и минимумов

Точки, где достигается наибольшее или наименьшее значение функции называются соответственно точками максимума или минимума функции.

Определение 1. Точка $$x_0$$ называется точкой максимума функции $$f(x)$$, если в некоторой окрестности точки $$x_0$$ выполняется неравенство $$f(x)\ge f(x_0)$$ (см. рис. 3.3).

Определение 2. Точка $$x_1$$ называется точкой минимума функции $$f(x)$$, если в некоторой окрестности точки $$x_1$$ выполняется неравенство $$f(x)\le f(x_1)$$ (см. рис. 3.3).

Значения функции в точках $$x_0$$ и $$x_1$$ называются соответственно максимумом и минимумом функции. Максимум и минимум функции объединяются общим названием экстремума функции.

Экстремум функции часто называют локальным экстремумом, подчеркивая тем самым, что понятие экстремума связано лишь с достаточно малой окрестностью точки $$x_0$$. Так что на одном промежутке функция может иметь несколько экстремумов, причем может случиться так, что минимум в одной точке больше максимума в другой.

Наличие максимума (или минимума) в отдельной точке промежутка $$X$$ вовсе не означает, что в этой точке функция $$f(x)$$ принимает наибольшее (наименьшее) значение на этом промежутке (или, как говорят имеет глобальный максимум (минимум)).

3.4.1.1 Теорема Ферма

Теорема Ферма. Если дифференцируемая на промежутке $$X$$ функция $$y = f(x)$$ достигает наибольшего или наименьшего значения в внутренней точке $$x_0$$, то тогда производная функции в этой точке равна нулю, т.е. $$f'(x_0) = 0$$.

(рис 3.3) Экстремумы функции

Пусть функция $$y = f(x)$$ дифференцируема на промежутке $$X$$ и в точке $$x_0\in X$$ принимает наименьшее значение (см. рис. 3.4).

Тогда $$f(x_0+\Delta x)\ge f(x_0)$$ если $$x_0+\Delta x\in X$$ и, следовательно $$\Delta y=f(x_0+\Delta x)-f(x_0)\ge 0$$ при достаточно малых $$\Delta x$$ и независимо от знака $$\Delta x$$.

Поэтому $$\frac{\Delta y}{\Delta x}\ge 0\; \mbox{при}\; \Delta x>0\; (\mbox{справа от}\; x_0);$$ $$\frac{\Delta y}{\Delta x}\le 0\; \mbox{при}\; \Delta x<0\; (\mbox{слева от}\; x_0).$$

(рис 3.4) Иллюстрация теоремы Ферма

Переходя к пределу справа и слева получим $$\lim_{\Delta x\to 0+}\frac{\Delta y}{\Delta x}\ge 0\; \mbox{и}\; \lim_{\Delta x\to 0-}\frac{\Delta y}{\Delta x}\le 0.$$ Так как функция дифференцируема на промежутке $$X$$, то пределы справа и слева равны $$\lim_{\Delta x\to 0+}\frac{\Delta y}{\Delta x}=\lim_{\Delta x\to 0-}\frac{\Delta y}{\Delta x}.$$ Отсюда $$f'(x_0) = 0$$.

Аналогичную последовательность рассуждений можно построить и для максимума.

Теорему Ферма часто называют необходимым условием экстремума дифференцируемой функции.

Геометрический смысл теоремы Ферма: в точке экстремума, достигаемого внутри промежутка $$X$$, касательная к графику функции параллельна оси абсцисс.

3.4.1.2 Необходимое условие экстремума

Если в точке $$x_0$$ дифференцируемая функция $$f(x)$$ имеет экстремум, то в некоторой окрестности этой точки выполняются условия теоремы Ферма, и следовательно, производная функции в этой точке равна нулю, т.е. $$f'(x_0) = 0$$. Но функция может иметь экстремум и в точках, в которых она не дифференцируема. Так, например, функция $$y = |x|$$ имеет экстремум (минимум) в точке $$x = 0$$, но не дифференцируема в ней. Функция $$y=\sqrt[3]{x^2}$$ также имеет в точке $$x = 0$$ минимум, а ее производная в этой точке бесконечна: $$y'=\displaystyle{\frac{2}{3\sqrt[3]{x}}},\ y'(0)=\infty$$.

Поэтому необходимое условие экстремума может быть сформулировано следующим образом.

Для того чтобы функция $$y = f(x)$$ имела экстремум в точке $$x_0$$, необходимо, чтобы ее производная в этой точке равнялась нулю ($$f'(x_0) = 0$$) или не существовала.

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

Пример. Найти критические точки функции и убедиться в наличии или отсутствии экстремума в этих точках: 1. $$y=x^2+1;$$ 2. $$y=x^3-1$$.

$$y' = 2x.\ y'(x) = 0$$ при $$x = 0$$. В точке $$x = 0$$ функция $$y = x^2 + 1$$ имеет минимум.

$$y' = 3x^2.\ y'(x) = 0$$ при $$x = 0$$. В точке $$x = 0$$ функция $$y = x^3-1$$ не имеет экстремума. Функция $$y = x^3-1$$ возрастает на всей числовой оси.

Итак, для нахождения экстремумов функции требуется дополнительное исследование критических точек.

Пример: Исследовать на наличие экстремума следующую функцию $$y(x) = x^3 - 3x^2 + 3x+2$$

Задаём исследуемую функцию

(%i1)	f(x):=x^3-3*x^2+3*x+2;

$$f\left( x\right) :={x}^{3}-3\,{x}^{2}+3\,x+2\leqno{(\%o1) }$$

Производную в форме функции определяем явно, используя функцию $$define$$

(%i2)	define(df(x),diff(f(x),x));

$$df\left( x\right) :=3\,{x}^{2}-6\,x+3\leqno{(\%o2) }$$

Решая уравнение $$df(x) = 0$$ (т.е. $$f'(x) = 0$$, находим критические точки

(%i3)	solve(df(x)=0,x);

$$[x=1]\leqno{(\%o3) }$$

В данном случае критическая точка одна — $$x = 1$$.

3.4.1.3 Первое достаточное условие экстремума

Теорема. Если при переходе через точку $$x_0$$ производная дифференцируемой функции $$y = f(x)$$ меняет свой знак с плюса на минус, то точка $$x_0$$ есть точка максимума функции $$y = f(x)$$, а если с минуса на плюс, то — точка минимума.

Пусть производная меняет знак с плюса на минус, т.е. в некотором интервале ($$a,x_0$$) производная положительна ($$f'(x) > 0$$), а в некотором интервале ($$x_0,b$$) — отрицательна ($$f'(x) < 0$$) (см. рис. 3.5). Тогда в соответствии с достаточным условием монотонности функция $$f(x)$$ возрастает на интервале ($$a,x_0$$) и убывает на интервале ($$x_0,b$$).

По определению возрастающей функции $$f(x_0)\ge f(x)$$ при всех $$x\in (a,x_0)$$, а по определению убывающей функции $$f(x)\le f(x_0)$$ при всех $$x\in (x_0,b)$$, т.е. $$f(x_0)\ge f(x)$$ при всех $$x\in (a,b)$$, следовательно, $$x_0$$ — точка максимума функции $$y = f(x)$$.

Аналогично рассматривается случай, когда производная меняет знак с минуса на плюс.

Отметим, что дифференцируемость функции в самой точке $$x_0$$ не использовалась при доказательстве теоремы. На самом деле она и не требуется — достаточно, чтобы функция была непрерывна в точке $$x_0$$.

Если изменение знака производной не происходит, то экстремума нет. Однако при работе с системами компьютерной математики удобнее второе достаточное условие экстремума.

(рис 3.5) Необходимое условие экстремума

3.4.1.4 Второе достаточное условие экстремума

Теорема. Если первая производная $$f'(x)$$ дважды дифференцируемой функции $$y = f(x)$$ равна нулю в некоторой точке $$x_0$$, а вторая производная в этой точке $$f''(x_0)$$ положительна, то $$x_0$$ есть точка максимума функции $$y = f(x)$$; если $$f''(x_0)$$ отрицательна, то $$x_0$$ — точка максимума.

Пусть $$f'(x_0)=0$$, а f$$f''(x_0)>0$$. Это значит, что $$f''(x)=(f'(x))'>0$$ также и в некоторой окрестности точки $$x_0$$, т.е. $$f'(x)$$ возрастает на некотором интервале $$(a,b)$$, содержащем точку $$x_0$$.

Но $$f'(x_0) = 0$$, следовательно, на интервале $$(a,x_0)\ f'(x)<0$$, а на интервале $$(x_0,b)\ f'(x)>0$$, т.е. $$f'(x)$$ при переходе через точку $$x_0$$ меняет знак с минуса на плюс, т.е. $$x_0$$ — точка минимума.

Аналогично рассматривается случай $$f'(x_0)=0$$ и $$f''(x_0)<0$$.

Продолжим исследование функции

$$y(x) = x^3 - 3x^2 + 3x + 2$$

Как установлено выше, имеется одна критическая точка: $$x = 1$$.

Задаёмся функцией $$d2f(x)$$

(%i4)	define(d2f(x),diff(df(x),x));

$$d2f\left( x\right) :=6\,x-6\leqno{(\%o4) }$$

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

(%i5)	map(d2f,%o3);

$$[6\,x-6=0]\leqno{(\%o5) }$$

В данном примере невозможно определить, является ли точка $$x = 1$$ экстремумом исследуемой функции, т.к. вторая производная в ней оказалась равной 0. Следует обратить внимание на способ вычисления — функция $$d2f(x)$$ применяется ко всем элементам списка, полученного при решении уравнения $$f'(x) = 0$$ (используется встроенная функция Maxima map).

Воспользуемся первым достаточным признаком наличия экстремума

(%i6)	df(0);

$$3\leqno{(\%o6) }$$

(%i7)	df(2);

$$3\leqno{(\%o7) }$$

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

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

(рис 3.6) Пример исследования функции

3.4.1.5 Схема исследования функции y = f(x) на экстремум

1. Найти производную $$y' = f'(x)$$.

2. Найти критические точки функции, в которых производная $$f'(x) = 0$$ или не существует.

3.1. Исследовать знак производной слева и справа от каждой критической точки и сделать вывод о наличии экстремумов функции.

 Или

3.2. Найти вторую производную $$f''(x)$$ и определить ее знак в каждой критической точке.

4. Найти экстремумы (экстремальные значения) функции.

Пример. Исследовать на экстремум функцию $$y = x(x - 1)^3$$.

1. y' = (x - 1)^3 + 3x(x - 1)^2 = (x - 1)^2(4x - 1).

2. Критические точки $$x_1 = 1$$ и $$x_2=\displaystyle{\frac{1}{4}}$$.

3. Изменение знака производной при переходе через точку $$x_1$$ не происходит, поэтому в этой точке нет экстремума.

$$y'' = 2(x - 1)(4x - 1) + 4(x - 1)^2 = 2[(x - 1)(6x - 3)]$$.

$$y''(x_2) > 0$$, поэтому в этой точке наблюдается минимум функции $$y = x(x - 1)^3$$.

4. $$y_{min}=y\left(\displaystyle{\frac{1}{4}}\right)=-\displaystyle{\frac{27}{256}}$$.

Выполним тот же расчёт при помощи Maxima

(%i13)	f(x):=x*(x-1)^3;

$$f\left( x\right) :=x\,{\left( x-1\right) }^{3}\leqno{(\%o13) }$$

(%i14)	define(df(x),diff(f(x),x));

$$df\left( x\right) :=3\,{\left( x-1\right) }^{2}\,x+{\left( x-1\right) }^{3}\leqno{(\%o14) }$$

(%i15)	solve(df(x)=0,x);

$$[x=\frac{1}{4},x=1]\leqno{(\%o15) }$$

(%i16)	define(d2f(x),diff(df(x),x));

$$d2f\left( x\right) :=6\,\left( x-1\right) \,x+6\,{\left( x-1\right) }^{2}\leqno{(\%o16) }$$

(%i17)	map(d2f,%o15);

$$[6\,\left( x-1\right) \,x+6\,{\left( x-1\right) }^{2}=\frac{9}{4},6\,\left( x-1\right) \,x+6\,{\left( x-1\right) }^{2}=0]\leqno{(\%o17) }$$

В точке $$x = 1$$ вторая производная равна 0, поэтому вычисляем значения первой производной слева и справа от $$x = 1$$ :

(%i18)	df(2);

$$7\leqno{(\%o18) }$$

(%i19)	df(1/3);

$$\frac{4}{27}\leqno{(\%o19) }$$

Производная в окрестности точки $$x = 1$$ не меняет знак, поэтому экстремум у исследуемой функции один — точка $$x=\frac{1}{4}$$. Так как $$d2f(\frac{1}{4}) > 0,\ x=\frac{1}{4}$$— точка минимума. Иллюстрация полученного результата — на рис. 3.7.

(рис 3.7) Пример исследования функции на экстремум

3.4.1.6 Нахождение наибольших и наименьших значений функции

Наибольшее или наименьшее значение функции на некотором отрезке может достигаться как в точках экстремума, так и в точках на концах отрезка.

Пусть функция $$y = f(x)$$ определена на некотором отрезке $$[a,b]$$.

Нахождение наибольших и наименьших значений функций происходит по следующей схеме.

1. Найти производную $$f'(x)$$.

2. Найти критические точки функции, в которых $$f'(x_0) = 0$$ или не существует.

3. Найти значения функции в критических точках и на концах отрезка и выбрать из них наибольшее $$f_{MAX}$$ и наименьшее $$f_{MIN}$$ значения. Это и будут наибольшее и наименьшее значение функции на исследуемом отрезке.

Пример. Найти наибольшее и наименьшее значения функции $$y =3x^2 - 6x$$ на отрезке [0, 3].

Аналитический расчёт:

1. $$y' = 6x - 6 ;y'' = 6$$.

2. $$x_0 = 1$$.

3. $$y(1) = -3 ; y(0) = 0 ; y(3) = 9$$.

В точке $$x = 1$$ наименьшее значение функции, а в точке $$x = 3$$ — наибольшее.

Расчёт с использованием Maxima:

Находим критические точки исследуемой функции

(%i29)	f(x):=3*x^2-6*x;

$$f\left( x\right) :=3\,{x}^{2}-6\,x\leqno{(\%o29) }$$

(%i30)	define(df(x),diff(f(x),x));

$$df\left( x\right) :=6\,x-6\leqno{(\%o30) }$$

(%i31)	solve(df(x)=0,x);

$$[x=1]\leqno{(\%o31) }$$

Результат расчёта — список, включающий один элемент ($$[x = 1]$$).

Создаём новый список, включающий граничные значений и критические точки:

(%i32)	L:[%o31[1],x=0,x=3];

$$[x=1,x=0,x=3]\leqno{(\%o32) }$$

Применяем функцию $$f(x)$$ к каждому элементу списка $$L$$:

(%i33)	map(f,L);

$$[3\,{x}^{2}-6\,x=-3,3\,{x}^{2}-6\,x=0,3\,{x}^{2}-6\,x=9]\leqno{(\%o33) }$$

Результат — наибольшие и наименьшие значения — находим в списке полученных значений.

3.4.2 Выпуклость функции

Определение. График функции $$y = f(x)$$ называется выпуклым в интервале $$(a,b)$$, если он расположен ниже касательной, проведенной в любой точке этого интервала (см. рис. 3.8а).

График функции $$y = f(x)$$ называется вогнутым в интервале $$(a,b)$$, если он расположен выше касательной, проведенной в любой точке этого интервала (см. рис. 3.8б).

(рис 3.8) Выпуклые и вогнутые функции.

3.4.2.1 Необходимые и достаточные условия выпуклости (вогнутости) функции

Для определения выпуклости (вогнутости) функции на некотором интервале можно использовать следующие теоремы.

Теорема 1. Пусть функция $$f(x)$$ определена и непрерывна на интервале $$X$$ и имеет конечную производную $$f'(x)$$. Для того, чтобы функция $$f(x)$$ была выпуклой (вогнутой) в $$X$$, необходимо и достаточно, чтобы ее производная $$f'(x)$$ убывала (возрастала) на этом интервале.

Теорема 2. Пусть функция $$f(x)$$ определена и непрерывна вместе со своей производной $$f'(x)$$ на $$X$$ и имеет внутри $$X$$ непрерывную вторую производную $$f''(x)$$. Для выпуклости (вогнутости) функции $$f(x)$$ в $$X$$ необходимо и достаточно, чтобы внутри $$X$$

$$f''(x)\le 0 ;f''(x)\ge 0.$$

Докажем теорему 2 для случая выпуклости функции $$f(x)$$.

Необходимость. Возьмем произвольную точку $$x_0\in X$$. Разложим функцию $$f(x)$$ около точки $$x_0$$ в ряд Тейлора

$$f(x)=f(x_0)+f'(x_0)(x-x_0)+r_1(x),$$ $$r_1(x)=\frac{(x-x_0)^2}{2}f''(x_0+\theta (x-x_0))\quad (0<\theta<1).$$

Уравнение касательной к кривой $$f(x)$$ в точке, имеющей абсциссу $$x_0$$:

$$Y(x)=f(x_0)+f'(x_0)(x-x_0).$$

Тогда превышение кривой $$f(x)$$ над касательной к ней в точке $$x_0$$ равно

$$f(x)-Y(x)=r_1(x).$$

Таким образом, остаток $$r_1(x)$$ равен величине превышения кривой $$f(x)$$ над касательной к ней в точке $$x_0$$. В силу непрерывности $$f''(x)$$, если $$f''(x_0) > 0$$, то и $$f''(x_0+\theta (x-x_0))>0$$ для $$x$$, принадлежащих достаточно малой окрестности точки $$x_0$$, а потому, очевидно, и $$r_1(x) > 0$$ для любого отличного от $$x_0$$ значения $$x$$, принадлежащего к указанной окрестности.

Значит, график функции $$f(x)$$ лежит выше касательной $$Y (x)$$ и кривая $$f(x)$$ выпукла в произвольной точке $$x_0\in X$$.

Достаточность. Пусть кривая $$f(x)$$ выпукла на промежутке $$X$$. Возьмем произвольную точку $$x_0\in X$$.

Аналогично предыдущему разложим функцию $$f(x)$$ около точки $$x_0$$ в ряд Тейлора

$$f(x)=f(x_0)+f'(x_0)(x-x_0)+r_1(x),$$ $$r_1(x)=\frac{(x-x_0)^2}{2}f''(x_0+\theta (x-x_0))\quad (0<\theta<1).$$

Превышение кривой $$f(x)$$ над касательной к ней в точке, имеющей абсциссу $$x_0$$, определяемой выражением $$Y(x)=f(x_0)+f'(x_0)(x-x_0)$$ равно

$$f(x)-Y(x)=r_1(x).$$

Так как превышение положительно для достаточно малой окрестности точки $$x_0$$, то положительна и вторая производная $$f''(x_0+\theta (x-x_0))$$. При стремлении $$x\to x_0$$ получаем, что для произвольной точки $$x_0$$ $$f''(x_0) > 0$$.

Пример. Исследовать на выпуклость (вогнутость) функцию $$y =x^2 - 16x + 32$$.

Ее производная $$y' = 2x - 16$$ возрастает на всей числовой оси, значит по теореме 1 функция вогнута на $$(-\infty,\infty)$$.

Ее вторая производная $$y'' = 2 > 0$$, поэтому по теореме 2 функция вогнута на $$(-\infty,\infty)$$.

3.4.2.2 Точки перегиба

Определение. Точкой перегиба графика непрерывной функции называется точка, разделяющая интервалы, в которых функция выпукла и вогнута.

Из этого определения следует, что точки перегиба — это точки точки экстремума первой производной. Отсюда вытекают следующие утверждения для необходимого и достаточного условий перегиба.

Теорема (необходимое условие перегиба). Для того чтобы точка $$x_0$$ являлась точкой перегиба дважды дифференцируемой функции $$y = f(x)$$, необходимо, чтобы ее вторая производная в этой точке равнялась нулю ($$f''(x_0) = 0$$) или не существовала.

Теорема (достаточное условие перегиба). Если вторая производная $$f''(x)$$ дважды дифференцируемой функции $$y = f(x)$$ при переходе через некоторую точку $$x_0$$ меняет знак, то $$x_0$$ есть точка перегиба.

Отметим, что в самой точке вторая производная $$f''(x_0)$$ может не существовать.

Геометрическая интерпретация точек перегиба иллюстрируется рис. 3.9

В окрестности точки $$x_1$$ функция выпукла и график ее лежит ниже касательной, проведенной в этой точке. В окрестности точки $$x_2$$ функция вогнута и график ее лежит выше касательной, проведенной в этой точке. В точке перегиба $$x_0$$ касательная разделяет график функции на области выпуклости и вогнутости.

3.4.2.3 Исследование функции на выпуклость и наличие точек перегиба

1. Найти вторую производную $$f''(x)$$.

2. Найти точки, в которых вторая производная $$f''(x) = 0$$ или не существует.

(рис 3.9) Точки перегиба.

3. Исследовать знак второй производной слева и справа от найденных точек и сделать вывод об интервалах выпуклости или вогнутости и наличии точек перегиба.

Пример. Исследовать функцию $$y(x) = 2x^3 - 6x^2 + 15$$ на выпуклость и наличие точек перегиба.

1. $$y' = 6x^2 - 12x; y'' = 12x - 12$$.

2. Вторая производная равна нулю при $$x_0 = 1$$.

3. Вторая производная $$y''(x)$$ меняет знак при $$x_0 = 1$$, значит точка $$x_0 = 1$$ — точка перегиба.

На интервале $$(-\infty,1)\ y''(x)<0$$, значит функция $$y(x)$$ выпукла на этом интервале.

На интервале $$(1,\infty)\ y''(x)>0$$, значит функция $$y(x)$$ вогнута на этом интервале.

3.4.2.4 Общая схема исследования функций и построения графика

При исследовании функции и построении ее графика рекомендуется использовать следующую схему:

Найти область определения функции.

Исследовать функцию на четность — нечетность. Напомним, что график четной функции симметричен относительно оси ординат, а график нечетной функции симметричен относительно начала координат.

Найти вертикальные асимптоты.

Исследовать поведение функции в бесконечности, найти горизонтальные или наклонные асимптоты.

Найти экстремумы и интервалы монотонности функции.

Найти интервалы выпуклости функции и точки перегиба.

Найти точки пересечения с осями координат.

Исследование функции проводится одновременно с построением ее графика.

Пример. Исследовать функцию $$y(x)=f(x)=\displaystyle{\frac{1+x^2}{1-x^2}}$$ и построить ее график.

1. Область определения функции — $$(-\infty,-1)\bigcup(-1,1)\bigcup(1,\infty)$$.

2. Исследуемая функция — четная $$y(x) = y(-x)$$, поэтому ее график симметричен относительно оси ординат.

3. Знаменатель функции обращается в ноль при $$x = \pm 1$$, поэтому график функции имеет вертикальные асимптоты $$x = -1$$ и $$x = 1$$.

Точки $$x = \pm 1$$ являются точками разрыва второго рода, так как пределы слева и справа в этих точках стремятся к $$\infty$$.

$$\lim_{x\to 1-0}y(x)=\lim_{x\to -1+0}y(x)=\infty; \lim_{x\to 1+0}y(x)=\lim_{x\to -1-0}y(x)=-\infty.$$

4. Поведение функции в бесконечности. $$\lim_{x\to\pm\infty}y(x)=-1,$$ поэтому график функции имеет горизонтальную асимптоту $$y = -1$$.

5. Экстремумы и интервалы монотонности. Находим первую производную

$$y'(x)=\frac{4x}{(1-x^2)}.$$

$$y'(x) < 0$$ при $$x\in(-\infty,-1)\bigcup(-1,0)$$, поэтому в этих интервалах функция $$y(x)$$ убывает.

$$y'(x) > 0$$ при $$x\in(0,1)\bigcup(1,\infty)$$, поэтому в этих интервалах функция $$y(x)$$ возрастает.

$$y'(x) = 0$$ при $$x = 0$$, поэтому точка $$x_0 = 0$$ является критической точкой.

Находим вторую производную

$$y''(x)=\frac{4(1+3x^2)}{(1-x^2)^3}.$$

Так как $$y''(0) > 0$$, то точка $$x_0 = 0$$ является точкой минимума функции $$y(x)$$.

6. Интервалы выпуклости и точки перегиба.

Функция $$y''(x) > 0$$ при $$x\in (-1, 1)$$, значит на этом интервале функция $$y(x)$$ вогнута.

Функция $$y''(x) < 0$$ при $$x\in(-\infty,-1)\bigcup(1,\infty)$$, значит на этих интервалах функция $$y(x)$$ выпукла.

Функция $$y''(x)$$ нигде не обращается в ноль, значит точек перегиба нет.

7. Точки пересечения с осями координат.

Уравнение $$f(0) = y$$, имеет решение $$y = 1$$, значит точка пересечения графика функции $$y(x)$$ с осью ординат (0, 1).

Уравнение $$f(x) = 0$$ не имеет решения, значит точек пересечения с осью абсцисс нет.

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

$$y(x)=\frac{1+x^2}{1-x^2}.$$

Схематически график функции $$y(x)=\frac{1+x^2}{1-x^2}$$ изображен на рис. 3.10.

(рис 3.10) График функции

3.4.2.5 Асимптоты графика функции

Определение. Асимптотой графика функции y = f(x) называется прямая, обладающая тем свойством, что расстояние от точки ($$x,f(x)$$) до этой прямой стремится к 0 при неограниченном удалении точки графика от начала координат.

Асимптоты бывают 3 видов: вертикальные (см. рис. 3.11а), горизонтальные (см. рис. 3.11б) и наклонные (см. рис. 3.11в).

Асимптоты находят, используя следующие теоремы:

Теорема 1. Пусть функция $$y = f(x)$$ определена в некоторой окрестности точки $$x_0$$ (исключая, возможно, саму эту точку) и хотя бы один из пределов функции при $$x \to x_0 \to 0$$ (слева) или $$x \to x_0 + 0$$ (справа) равен бесконечности. Тогда прямая является $$x = x_0$$ вертикальной асимптотой графика функции $$y = f(x)$$.

Вертикальные асимптоты $$x = x_0$$ следует искать в точках разрыва функции $$y = f(x)$$.

Теорема 2. Пусть функция $$y = f(x)$$ определена при достаточно больших $$x$$ и существует конечный предел функции $$\lim_{x\to\mp\infty}f(x)=b.$$ Тогда прямая $$y = b$$ есть горизонтальная асимптота графика функции $$y = f(x)$$.

(рис 3.11) Асимптоты

Теорема 3. Пусть функция $$y = f(x)$$ определена при достаточно больших $$x$$ и существуют конечные пределы $$\lim_{x\to\mp\infty}\frac{f(x)}{x}=k$$ и $$\lim_{x\to\mp\infty}[f(x)-kx]=b.$$ Тогда прямая $$y = kx + b$$ является наклонной асимптотой графика функции $$y = f(x)$$.

Пример. Найти асимптоты графика дробно-рациональной функции

$$y(x)=\frac{ax+b}{cx+d}; c\ne 0; ad-bc\ne 0.$$

Если $$c = 0$$, то дробно-рациональная функция становится линейной

$$y(x)=\frac{a}{d}x+\frac{b}{d}.$$

Особая точка $$x = -d/c$$. Найдём предел $$\displaystyle{\lim_{x\to -d/c}f(x)}$$.

Перепишем дробно-рациональную функцию в виде: $$y(x)=\frac{ax+b}{c(x+d/c)}$$ Так как $$ad - bc \ne 0$$ то при $$x \to d/c$$ числитель дробно-рациональной функции не стремится к нулю. Поэтому прямая $$x = -d/c$$ — асимптота графика дробно-рациональной функции.

Найдём предел $$\displaystyle{\lim_{x\to \pm\infty}f(x)}$$. $$\lim_{x\to \pm\infty}\frac{ax+b}{cx+d}= \lim_{x\to \pm\infty}\frac{a+b/x}{c+d/c}=\frac{a}{c}$$ $$y = a/c$$ — является горизонтальной асимптотой дробно-рациональной функции.

Пример. Найти асимптоты кривой $$y(x)=\displaystyle{\frac{x^3}{x^2+1}}$$. $$\lim_{x\to\pm\infty}\frac{f(x)}{x}= \lim_{x\to\pm\infty}\frac{x^2}{x^2+1}=1.$$ Поэтому $$k = 1$$.

Теперь ищем $$b$$.

$$b=\lim_{x\to\pm\infty}\left[\frac{x^3}{x^2+1}-x\right]= \lim_{x\to\pm\infty}\left(\frac{-x}{x^2+1}\right)$$

Функция $$y(x)=\displaystyle{\frac{x^3}{x^2+1}}$$ имеет наклонную асимптоту $$y = x$$.

3.4.2.6 Свойства функций, непрерывных на отрезке. Теоремы Вейерштрасса

Если функция $$y = f(x)$$ непрерывна на отрезке $$[a,b]$$, то она ограничена на этом отрезке, т.е. существуют такие постоянные и конечные числа $$m$$ и $$M$$, что $$m\le f(x)\le M \quad \text{при}\quad a\le x\le b$$ (см. рис. 3.12а).

Если функция $$y = f(x)$$ непрерывна на отрезке [$$a,b$$], то она достигает на этом отрезке наибольшего значения $$M$$ и наименьшего значения m (см. рис. 3.12б).

Если функция $$y = f(x)$$ непрерывна на отрезке [$$a,b$$], и значения её на концах отрезка $$f(a)$$ и $$f(b)$$ имеют противоположные знаки, то внутри отрезка найдётся точка $$\xi\in(a,b)$$, такая, что $$f(\xi)=0$$ (см. рис. 3.12в).

(рис 3.12) Иллюстрации к теоремам Вейерштрасса

3.4.3 Дифференцирование функций нескольких переменных

Для определения набора частных производных функции нескольких переменным (компонентов градиента) используется функция $$gradef$$ в формате $$gradef(f(x_1,...,x_n),g_1,...,g_m)$$ или $$gradef(a,x,expr)$$

Выражение $$gradef(f(x_1,...,x_n),g_1,...,g_m)$$ определяет $$g_1,g_2,... g_n$$ как частные производные функции $$f(x_1,x_2,... ,x_n)$$ по переменным $$x_1,x_2,... ,x_n$$ соответственно.

Зависимости между переменными можно явно указать при помощи функции $$depends$$, которая позволяет декларировать, что переменная зависит от одной или нескольких других переменных. Например, если зависимость $$f$$ и $$x$$ отсутствует, выражение $$diff(f,x)$$ возвращает 0. Если декларировать её при помощи $$depends(f,x)$$, выражение $$diff(f,x)$$ возвращает символьную производную.

Пример:

(%i1)	depends(y,x);

$$[y\left( x\right) ]\leqno{(\%o1) }$$

(%i2)	gradef(f(x,y),x^2,g(x,y));

$$f\left( x,y\right) \leqno{(\%o2) }$$

(%i3)	diff(f(x,y),x);

$$g\left( x,y\right) \,\left( \frac{d}{d\,x}\,y\right) +{x}^{2}\leqno{(\%o3) }$$

(%i4)	diff(f(x,y),y);

$$g\left( x,y\right) \leqno{(\%o4) }$$

Вторая форма обращения к $$gradef$$ фактически устанавливает зависимость $$a$$ от $$x$$. При помощи $$gradef$$ можно определить производные некоторой функции, даже если она сама неизвестна, посредством $$diff$$ определить производные высших порядков.

Для прямых вычислений, связанных с операциями векторного анализа, необходимо загрузить пакет vect. Кроме того, применения операторов $$div, curl, grad, laplasian$$ к некоторому выражению используется функция $$express$$.

Пример: Вычисление градиента функции трех переменных

(%i2)	grad (x^2 + 2*y^2 + 3*z^2);

$$grad\left( 3\,{z}^{2}+2\,{y}^{2}+{x}^{2}\right) \leqno{(\%o2) }$$

(%i3)	express(%);

$$[\frac{d}{d\,x}\,\left( 3\,{z}^{2}+2\,{y}^{2}+{x}^{2}\right) ,\frac{d}{d\,y}\,\left( 3\,{z}^{2}+2\,{y}^{2}+{x}^{2}\right) ,\frac{d}{d\,z}\,\left( 3\,{z}^{2}+2\,{y}^{2}+{x}^{2}\right) ]\leqno{(\%o3) }$$

(%i4)	ev(%,diff);

$$[2\,x,4\,y,6\,z]\leqno{(\%o4) }$$

Вычисление дивергенции

(%i5)	div([x^2,2*y^2,3*z^2]);

$$div\left( [{x}^{2},2\,{y}^{2},3\,{z}^{2}]\right) \leqno{(\%o5) }$$

(%i6)	express(%);

$$\frac{d}{d\,z}\,\left( 3\,{z}^{2}\right) +\frac{d}{d\,y}\,\left( 2\,{y}^{2}\right) +\frac{d}{d\,x}\,{x}^{2}\leqno{(\%o6) }$$

(%i7)	ev(%,diff);

$$6\,z+4\,y+2\,x\leqno{(\%o7) }$$

Вычисление вихря:

(%i8)	curl([x^2,2*y^2,3*z^2]);

$$curl\left( [{x}^{2},2\,{y}^{2},3\,{z}^{2}]\right) \leqno{(\%o8) }$$

(%i9)	express(%);

$$[\frac{d}{d\,y}\,\left( 3\,{z}^{2}\right) -\frac{d}{d\,z}\,\left( 2\,{y}^{2}\right) ,\frac{d}{d\,z}\,{x}^{2}-\frac{d}{d\,x}\,\left( 3\,{z}^{2}\right) ,\frac{d}{d\,x}\,\left( 2\,{y}^{2}\right) -\frac{d}{d\,y}\,{x}^{2}]\leqno{(\%o9) }$$

(%i10)	ev(%,diff);

$$[0,0,0]\leqno{(\%o10) }$$

Вычисление оператора Лапласа:

(%i13)	laplacian(x^2+2*y^2+3*z^2);

$$laplacian\left( 3\,{z}^{2}+2\,{y}^{2}+{x}^{2}\right) \leqno{(\%o13) }$$

(%i14)	express(%);

$$\frac{{d}^{2}}{d\,{z}^{2}}\,\left( 3\,{z}^{2}+2\,{y}^{2}+{x}^{2}\right) +\frac{{d}^{2}}{d\,{y}^{2}}\,\left( 3\,{z}^{2}+2\,{y}^{2}+{x}^{2}\right) +\frac{{d}^{2}}{d\,{x}^{2}}\,\left( 3\,{z}^{2}+2\,{y}^{2}+{x}^{2}\right) \leqno{(\%o14) }$$

(%i15)	ev(%,diff);

$$12\leqno{(\%o15) }$$

Рассмотрим пример исследования функции нескольких переменных: исследовать на экстремум функцию $$f(x,y) = {y}^{2}-4\,y+{x}^{3}-\frac{9\,{x}^{2}}{2}+6\,x-12$$

Загружаем пакет vect

(%i1)	load("vect")$

Определяем исследуемое выражение и вычисляем его градиент:

(%i2)	f:x^3-9/2*x^2+6*x+y^2-4*y-12;

$${y}^{2}-4\,y+{x}^{3}-\frac{9\,{x}^{2}}{2}+6\,x-12\leqno{(\%o2) }$$

(%i3)	grad(f);

$$grad\left( {y}^{2}-4\,y+{x}^{3}-\frac{9\,{x}^{2}}{2}+6\,x-12\right) \leqno{(\%o3) }$$

(%i4)	express(%);

$$[\frac{d}{dx}\left( {y}^{2}-4y+{x}^{3}-\frac{9{x}^{2}}{2}+6x-12\right),\leqno{(\%o4) }$$ $$\frac{d}{d\,y}\left( {y}^{2}-4y+{x}^{3}-\frac{9{x}^{2}}{2}+6x-12\right),$$ $$\frac{d}{dz}\,\left( {y}^{2}-4y+{x}^{3}-\frac{9{x}^{2}}{2}+6x-12\right) ]$$

(%i5)	ev(%,diff);

$$[3\,{x}^{2}-9\,x+6,2\,y-4,0]\leqno{(\%o5) }$$

Выделяем из полученного списка частные производные и решаем систему $$f_x(x,y) = 0;\,f_y(x,y) = 0$$

(%i6)	dfdx:%o5[1];

$$3\,{x}^{2}-9\,x+6\leqno{(\%o6) }$$

(%i7)	dfdy:%o5[2];

$$2\,y-4\leqno{(\%o7) }$$

(%i8)	solve([dfdx=0,dfdy=0],[x,y]);

$$[[x=1,y=2],[x=2,y=2]]\leqno{(\%o8) }$$

В результате решения находим две критические точки $$M_1(1, 2)$$ и $$M_2(2, 2)$$ Для проверки, достигается ли в критических точках экстремум, используем достаточное условие экстремума:

(%i9)	A:diff(dfdx,x);

$$6\,x-9\leqno{(\%o9) }$$

(%i10)	C:diff(dfdy,y);

$$2\leqno{(\%o10) }$$

(%i11)	B:diff(dfdx,y);

$$0\leqno{(\%o11) }$$

(%i12)	A*C-B^2;

$$2\,\left( 6\,x-9\right) \leqno{(\%o12) }$$

Так как $$A * C - B^2 > 0$$ только в точке $$M_2(2, 2)$$, то исследуемая функция имеет единственный экстремум. Учитывая, что в точке $$M_2(2, 2) A > 0$$, точка $$M_2$$ — точка минимума. Результат иллюстрируем графически рис. 3.13).

(рис 3.13) Поиск экстремума функции нескольких переменных

3.5 Аналитическое и численное интегрирование

3.5.1 Основные команды

Неопределённый интеграл $$\int f(x)dx$$ вычисляется с помощью команды $$integrate(f, x)$$, где $$f$$ — подинтегральная функция, $$x$$ — переменная интегрирования.

Для вычисления определённого интеграла $$\int _{a}^{b}f(x)dx$$ в команде $$integrate$$ добавляются пределы интегрирования, например, $$integrate((1+cos(x))^2, x, 0, %pi);$$

$$\[\int _{0}^{\pi }(1+\cos (x))^{2} dx=\frac{3}{2} \pi \]$$

Несобственные интегралы с бесконечными пределами интегрирования вычисляются, если в параметрах команды $$integrate$$ указывать, например, $$x, 0, inf$$.

Численное интегрирование выполняется функцией $$romberg$$ или при помощи функций пакета quadpack.

3.5.2 Интегралы, зависящие от параметра. Ограничения для параметров

Если требуется вычислить интеграл, зависящий от параметра, то его значение может зависеть от знака этого параметра или каких- либо других ограничений. Рассмотрим в качестве примера интеграл $$\int _{0}^{+\infty }e^{-ax} dx$$, который, как известно из математического анализа, сходится при $$а>0$$ и расходится при $$а<0$$. Если вычислить его сразу, то получится:

(%i1)	integrate(exp(-a*x),x,0,inf);
Is a positive, negative, or zero? p;

$$(\%o1)\ \frac{1}{a}$$

Результат аналитического интегрирования

$$\[\int _{0}^{+\infty }e^{(-ax)} dx=\mathop{\lim }\limits_{x\to \infty } -\frac{e^{(-ax)} -1}{a} .\]$$

Для получения явного аналитического результата вычислений следует сделать какие-либо предположения о значении параметров, то есть наложить на них ограничения. Это можно сделать при помощи команды $$assume(expr1)$$, где $$expr1$$ — неравенство.

Описание наложенных ограничений параметра a можно вызвать командой $$properties(a)$$.

(%i1)	assume (a > 1)$ integrate (x**a/(x+1)**(5/2), x, 0, inf);
Is (2a+2)/5 an integer? no;
Is 2a-3 positive, negative, or zero? neg;

$$(\%o2)\ \beta\left( a+1,\frac{3}{2}-a\right)$$

(%i3)	properties(a);

$$(\%o3)\ [\text{database info, }\ a>1]$$

Вернёмся к вычислению интеграла с параметром $$\int _{0}^{+\infty }e^{-ax} dx$$, которое следует производить в таком порядке:

(%i1)	assume(a>0); integrate(exp(-a*x),x,0,inf);

$$(\%o1)\ [a>0]\quad (\%o2)\ \frac{1}{a}$$

Отменить принятые ограничения на значения параметров можно, используя функцию $$forget$$.

Пример:

(%i1)	assume(n+1>0); integrate((a+b)*x^(n+1),x);

$$(\%o1)\ [n>-1]\quad (\%o2)\ \frac{\left( b+a\right) \,{x}^{n+2}}{n+2}$$

Отмена ограничения влечёт за собой вопрос о значениях параметров подинтегральной функции:

(%i3)	forget(n+1>0); integrate((a+b)*x^(n+1),x);

$$(\%o3)\quad [n>-1]$$

Is n + 2 zero or nonzero? zero;

$$(\%o4)\quad \left( b+a\right) \,log\left( x\right)$$

Результат, который получен, совершенно другой!

3.5.3 Основные приёмы интегрирования

В Maxima имеется функция, предназначенных для выполнения расчётов шаг за шагом, осуществляющая замену переменной $$changevar$$.

Формулу интегрирования по частям: $$\[\int u(x)v'(x)dx=u(x)v(x)-\int u'(x)v(x)dx \]$$ придётся применять вручную. В Maxima (в отличие от, например, Maple), функция интегрирования по частям не выделена явно, хотя в отдельных случаях этот способ используется $$integrate$$.

Для вычисления первообразных дифференциальных выражений используется пакет antid (основные функции пакета — $$antidiff$$ и $$antid$$). Функция $$antidiff$$ выполняет интегрирование выражений с произвольными функциями (в том числе неопределёнными), перед её первым вызовом следует загрузить пакет ($$antid$$ отличается от неё форматом выводимого результата).

Пример:

(%i1)	load("antid");
(%i2)	expr: exp(z(x))*diff(z(x),x)*sin(x);

$$(\%o2)\ {e}^{z\left( x\right) }\,sin\left( x\right) \,\left( \frac{d}{d\,x}\,z\left( x\right) \right)$$

(%i3)	a1: antid (expr, x, z(x));

$$(\%o3)\ [{e}^{z\left( x\right) }\,sin\left( x\right) ,-{e}^{z\left( x\right) }\,cos\left( x\right) ]$$

При помощи пакета antid можно выполнить формальное интегрирование по частям, например:

(%i1)	expr:u(x)*diff(v(x),x);

$$(\%o1)\ \mathrm{u}\left( x\right) \,\left( \frac{d}{d\,x}\,\mathrm{v}\left( x\right) \right)$$

(%i2)	a:antid(expr,x,v(x));

$$(\%o2)\ \mathrm{antid}\left( \mathrm{u}\left( x\right) \,\left( \frac{d}{d\,x}\,\mathrm{v}\left( x\right) \right) ,x,\mathrm{v}\left( x\right) \right)$$

(%i3)	b:antidiff(expr,x,v(x));

$$(\%o3)\ \mathrm{antidiff}\left( \mathrm{u}\left( x\right) \,\left( \frac{d}{d\,x}\,\mathrm{v}\left( x\right) \right) ,x,\mathrm{v}\left( x\right) \right)$$

Если в интеграле требуется сделать замену переменных, используется функция $$changevar$$.

Синтаксис вызова этой функции: $$changevar(expr,f(x,y),y,x)$$.

Функция осуществляет замену переменной в соответствии с уравнением $$f(x,y) = 0$$ во всех интегралах, встречающихся в выражении expr (предполагается, что $$y$$ — новая переменная, $$x$$ — исходная). При использовании совместно с $$changevar$$ часто используется отложенное вычисление интеграла (одинарная кавычка перед функцией $$integrate$$).

Пример:

(%i5)	assume(a > 0)$ 'integrate (%e**sqrt(a*y), y, 0, 4);

$$\int_{0}^{4}{e}^{\sqrt{a}\,\sqrt{y}}dy\leqno{(\%o6) }$$

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

(%i7)	changevar (%, y-z^2/a, z, y);

$$-\frac{2\,\int_{-2\,\sqrt{a}}^{0}z\,{e}^{\left| z\right| }dz}{a}\leqno{(\%o7) }$$

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

(%i8)	ev(%,nouns);

$$-\frac{2\,\left( -2\,\sqrt{a}\,{e}^{2\,\sqrt{a}}+{e}^{2\,\sqrt{a}}-1\right) }{a}\leqno{(\%o8) }$$

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

Пример:

(%i10)	expand ((x-4) * (x^3+2*x+1));

$${x}^{4}-4\,{x}^{3}+2\,{x}^{2}-7\,x-4\leqno{(\%o10) }$$

(%i11)	integrate (1/%, x);

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

$$\frac{log\left( x-4\right) }{73}-\frac{\int \frac{{x}^{2}+4\,x+18}{{x}^{3}+2\,x+1}dx}{73}\leqno{(\%o11) }$$

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

Рассмотрим предыдущий пример, выполнив предварительно факторизацию знаменателя:

(%i1)	f:expand ((x-4) * (x^3+2*x+1));

$${x}^{4}-4\,{x}^{3}+2\,{x}^{2}-7\,x-4\leqno{(\%o1) }$$

(%i2)	polyfactor:true$ ffact:allroots(f);

$$1.0 ( x-3.999999999999997)( x+0.4533976515164)\leqno{(\%o3)}$$ $$( {x}^{2}-0.45339765151641\,x+2.205569430400593)$$

(%i4)	float(integrate(1/ffact,x));

Полученный результат всё равно трудно назвать однозначно приемлемым, т.к. он включает одновременно очень большие и очень малые величины. Причина в том, что корни знаменателя представлялись рациональными числами. Для того, чтобы получить компактный результат, желательно для коэффициентов вида $$r=\frac{m}{n}$$ уменьшить $$m$$ и $$n$$.

Интегралы от тригонометрических и логарифмических функций Maxima вычисляет довольно успешно. Рассмотрим несколько примеров.

(%i1)	integrate(sin(x)*sin(2*x)*sin(3*x),x);

$$\frac{cos\left( 6\,x\right) }{24}-\frac{cos\left( 4\,x\right) }{16}-\frac{cos\left( 2\,x\right) }{8}\leqno{(\%o1) }$$

(%i2)	integrate(1/cos(x)^3,x);

$$\frac{log\left( sin\left( x\right) +1\right) }{4}-\frac{log\left( sin\left( x\right) -1\right) }{4}-\frac{sin\left( x\right) }{2\,{sin\left( x\right) }^{2}-2}\leqno{(\%o2) }$$

(%i3)	integrate(x^3*log(x),x);

$$\frac{{x}^{4}\,log\left( x\right) }{4}-\frac{{x}^{4}}{16}\leqno{(\%o3) }$$

3.5.4 Преобразование Лапласа

Прямое и обратное преобразование Лапласа вычисляются посредством функций $$laplace$$ и $$ilt$$ соответственно.

Синтаксис обращения к функции $$laplace: laplace(expr,t,s)$$.

Функция вычисляет преобразование Лапласа выражения $$expr$$ по отношению к переменной $$t$$. Образ выражения $$expr$$ будет включать переменную $$s$$.

Функция $$laplace$$ распознаёт в выражении expr функции $$delta, exp, log, sin, cos, sinh, cosh$$, и $$erf$$, а также производные, интегралы, суммы и обратное преобразование Лапласа ($$ilt$$). При наличии других функций вычисление преобразования может и не удаться.

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

(%i1)	laplace(c,t,s);

$$\frac{c}{s}\leqno{(\%o1) }$$

(%i2)	laplace(erf(t),t,s);

$$\frac{{e}^{\frac{{s}^{2}}{4}}\,\left( 1-erf\left( \frac{s}{2}\right) \right) }{s}\leqno{(\%o2) }$$

(%i3)	laplace(sin(t)*exp(-a*t),t,s);

$$\frac{1}{{s}^{2}+2\,a\,s+{a}^{2}+1}\leqno{(\%o3) }$$

Функция $$ilt(expr,t,s)$$ вычисляет обратное преобразование Лапласа относительно переменной $$t$$ с параметром $$s$$.

Пример:

(%i1)	laplace(c,t,s);

$$\frac{c}{s}\leqno{(\%o1) }$$

(%i2)	ilt(%,s,t);

$$c\leqno{(\%o2) }$$

(%i3)	laplace(sin(2*t)*exp(-4*t),t,s);

$$\frac{2}{{s}^{2}+8\,s+20}\leqno{(\%o3) }$$

(%i4)	ilt(%,s,t);

$${e}^{-4\,t}\,sin\left( 2\,t\right) \leqno{(\%o4) }$$

3.6 Методы теории приближения в численном анализе

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

Этим объясняется важность всевозможных аналитических и численных приложений полиномиальных приближений для аппроксимации и вычисления функции. Замена функций на их степенные разложения и полиномиальные приближения помогает изучению пределов, анализу сходимости и расходимости рядов и интегралов, приближённому вычислению интегралов и решению дифференциальных уравнений. Степенные ряды и разложения по многочленам Чебышёва широко используются при вычислении значений функции с заданной степенью точности. Они являются эффективным вычислительным средством при решении широкого круга научно-технических задач.

3.6.1 Приближённое вычисление математических функций

Пусть функция $$f(x)$$ задана на интервале $$(x_0 - R,x_0 + R)$$ и нам требуется вычислить значение функции $$f(x)$$ при $$x = x_1 \in (x_0-R,x_0+R)$$ с заданной точностью $$\epsilon > 0$$.

Предположив, что функция $$f(x)$$ в интервале $$x\in (x_0 - R,x_0 + R)$$ раскладывается в степенной ряд $$f(x)=\sum^{\infty }_{i=0}u_{i}(x)= \sum^{\infty }_{i=0}a_{i}(x-x_{0})^{i}=a_{0}+ a_{1}(x-x_{0})+a_{2}(x-x_{0})^{2}+\dots +a_{n}(x-x_{0})^{n}+\dots,$$ мы получим, что точное значение $$f(x_1)$$ равно сумме этого ряда при $$x = x_1$$ $$f(x_{1})=\sum^{\infty }_{i=0}a_{i}(x_{1}-x_{0})^{i}=a_{0}+ a_{1}(x_{1}-x_{0})+a_{2}(x_{1}-x_{0})^{2}+\dots+ a_{n}(x_{1}-x_{0})^{n}+\dots,$$ а приближённое — частичной сумме $$S_n(x_1)$$ $$f(x_{1})\approx S_{n}(x_{1})=\sum^{n}_{i=0}a_{i}(x_{1}- x_{0})^{i}=a_{0}+a_{1}(x_{1}-x_{0})+a_{2}(x_{1}-x_{0})^{2}+ \dots+a_{n}(x_{1}-x_{0})^{n}.$$

Для погрешности приближения мы имеем выражение в виде остатка ряда $$f(x_{1})-S_{n}(x_{1})=r_{n}(x_{1}),$$ где $$r_{n}(x_{1})=\sum^{\infty}_{i=1}x_{1}^{n+i}= a_{n+1}x_{1}^{n+1}+a_{n+2}x_{1}^{n+2}+\dots$$

Для знакопеременных рядов с последовательно убывающими членами $$\vert r_{n}(x)\vert=\vert \sum^{\infty}_{i=1}u_{n+i}(x_{1})\vert <\vert u_{n+1}(x_{1})\vert .$$

Точность аппроксимации, как правило, возрастает с ростом степени приближающего степенного разложения и тем выше, чем точка $$x$$ ближе к точке $$x_0$$. Для равномерной аппроксимации на интервале наиболее удобными оказываются разложения по многочленам Чебышёва.

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

Ряд Тейлора для функции $$f(x)$$ — это степенной ряд вида $$\sum^{\infty }_{k=0}\frac{f^{(k)}(x_{0})}{k!}(x-x_{0})^{k},$$ где числовая функция $$f$$ предполагается определённой в некоторой окрестности точки $$x_0$$ и имеющей в этой точке производные всех порядков.

Многочленами Тейлора для функции $$f(x)$$, порядка $$n$$ соответственно, называются частные суммы ряда Тейлора $$\sum^{n}_{k=0}\frac{f^{(k)}(x_{0})}{k!}(x-x_{0})^{k}.$$

Если мы распишем эту формулу, то получим следующее выражение $$f(x_{0})+\frac{f^{\prime}(x_{0})}{1!}(x-x_{0})+ \frac{f^{\prime\prime}(x_{0})}{2!}(x-x_{0})^{2}+\dots+ \frac{f^{(n)}(x_{0})}{n!}(x-x_{0})^{n}.$$

Формула Тейлора для функции $$f(x)$$ — это представление функции в виде суммы её многочлена Тейлора степени $$n(n = 0, 1, 2,...)$$ и остаточного члена. Другими словами это называют разложением функции $$f(x)$$ по формуле Тейлора в окрестности точки $$x_0$$. Если действительная функция $$f$$ одного переменного имеет $$n$$ производных в точке $$x_0$$, то её формула Тейлора имеет вид $$f(x)=P_{n}(x)+r_{n}(x),$$ где $$P_{n}(x)=\sum^{n}_{k=0}\frac{f^{(k)}(x_{0})}{k!}(x-x_{0})^{k}$$ — многочлен Тейлора степени $$n$$, а остаточный член может быть записан в форме Пеано $$r_{n}(x)=o((x-x_{0})^{n}), x\rightarrow x_{0}.$$

Получаем, что $$P_{n}(x)=f(x_{0})+\frac{f^{\prime}(x_{0})}{1!}(x-x_{0})+ \frac{f^{\prime\prime}(x_{0})}{2!}(x-x_{0})^{2}+\dots+ \frac{f^{(n)}(x_{0})}{n!}(x-x_{0})^{n}.$$

Если функция $$f$$ дифференцируема $$n+1$$ раз в некоторой окрестности точки $$x_{0}, (x_{0}-\delta , x_{0}+\delta ),\delta >0$$, то остаточный член в этой окрестности может быть записан в форме Лагранжа $$r_{n}(x)=\frac{f^{(n+1)}(x_{0}+\theta (x-x_{0}))}{(n+1)!}(x-x_{0})^{(n+1)} ,$$ $$0< \theta <1, x\in (x_{0}-\delta ,x_{0}+\delta ).$$

Заметим, что при $$n = 1$$ выражение для $$P_1(x) = f(x_0)+ f'(x_0)(x-x_0)$$ совпадает с формулой Лагранжа конечных приращений для функции $$f(x)$$.

Формула Тейлора для многочленов. Пусть имеется произвольный многочлен $$f(x)=a_{0}x^{n}+a_{1}x^{n-1}+\dots+a_{n}$$.. Тогда при любых $$x$$ и $$h$$ имеет место следующая формула:

$$f(x+h)=a_{0}(x+h)^{n}+a_{1}(x+h)^{n-1}+\dots+a_{n}=\\ =f(x)+f^{\prime}(x)h+\frac{f^{\prime\prime}(x)}{2!}h^{2}+\dots+ \frac{f^{(k)}(x)}{k!}h^{k}+\dots+\frac{f^{(n)}(x)}{n!}h^{n}.$$

Рядом Маклорена для функции $$f(x)$$ называется её ряд Тейлора в точке 0 начала координат, то есть таким образом это степенной ряд вида $$f(x)=\sum^{\infty }_{k=0}\frac{f^{(k)}(0)}{k!}x^{k} .$$

Таким образом формула Маклорена является частным случаем формулы Тейлора. Предположим, что функция $$f(x)$$ имеет n производных в точке $$x = 0$$ . Тогда в некоторой окрестности этой точки $$(-\delta , \delta ), \delta >0$$ , функцию $$f(x)$$ можно представить в виде $$f(x)=\sum^{n}_{k=0}\frac{f^{(k)}(0)}{k!}x^{k}+r_{n}(x),$$ $$x\in (-\delta , \delta ),$$ где $$r_n(x)$$ — остаточный член $$n$$-ого порядка в форме Пеано.

Приведём разложения по формуле Маклорена для основных элементарных математических функций:

$$e^{x}=1+x+\frac{x^{2}}{2!}+\frac{x^{3}}{3!}+\dots+\frac{x^{n}}{n!}+o(x^{n}),\\ sinx=x-\frac{x^{3}}{3!}+\frac{x^{5}}{5!}+\dots+ (-1)^{n-1}\frac{x^{2n-1}}{(2n-1)!}+o(x^{2n}),\\ cosx=1-\frac{x^{2}}{2!}+\frac{x^{4}}{4!}+\dots+ (-1)^{n}\frac{x^{2n}}{(2n)!}+o(x^{2n+1}),\\ (1+x)^{\alpha }=1+\alpha x+\frac{\alpha (\alpha -1)}{2!}x^{2}+ \dots+\frac{\alpha (\alpha -1)\dots(\alpha-n+1)}{n!}x^{n}+o(x^{n}),\\ ln(1+x)=x-\frac{x^{2}}{2}+\frac{x^{3}}{3}+\dots+ (-1)^{n-1}\frac{x^{n}}{n}+o(x^{n}).$$

В Maxima существует специальная команда, позволяющая вычислять ряды и многочлены Тейлора: $$taylor(expr,x,a,n)$$. Здесь $$expr$$ — разлагаемое в ряд выражение, $$a$$ — значение $$x$$, в окрестности которого выполняется разложение (по степеням $$x-a$$), $$n$$ — параметр, указывающий на порядок разложения и представленный целым положительным числом. Если $$a$$ указывается просто в виде имени переменной, то производится вычисление ряда и многочлена Маклорена.

Пример: Найти многочлен Тейлора 9-ой степени экспоненциальной функции $$e^x$$ в начале координат.

(%i29)	taylor(exp(x),x,0,9);

$$1+x+\frac{{x}^{2}}{2}+\frac{{x}^{3}}{6}+\frac{{x}^{4}}{24}+\frac{{x}^{5}}{120}+\frac{{x}^{6}}{720}+\frac{{x}^{7}}{5040}+\frac{{x}^{8}}{40320}+\frac{{x}^{9}}{362880}+\dots\leqno{(\%o29) }$$

Многочлены Тейлора дают наиболее точную аппроксимацию приближаемой функции вблизи точки $$x_0$$. По мере удаления от точки $$x_0$$ погрешность возрастает. Для приближения приходится использовать многочлены Тейлора более высокой степени, но иногда и они не помогают в связи с накоплением вычислительной погрешности.

Интересно проследить этот процесс графически. Пакет Maxima предоставляет такую возможность с помощью команды $$plot$$.

Пример: Найти число $$e$$ с точностью до 0.001. Положим $$x =1$$. Тогда чтобы вычислить значение $$e$$, необходимо выполнить серию команд:

Строим разложение функции $$e^x$$ в ряд Тейлора (до 8 порядка включительно)

(%i1)	t:taylor(exp(x),x,0,8);

$$1+x+\frac{{x}^{2}}{2}+\frac{{x}^{3}}{6}+\frac{{x}^{4}}{24}+\frac{{x}^{5}}{120}+\frac{{x}^{6}}{720}+\frac{{x}^{7}}{5040}+\frac{{x}^{8}}{40320}+\dots\leqno{(\%o1) }$$

Вычисляем частичную сумму ряда при $$x = 1$$:

(%i2)	ev(t,x=1);

$$\frac{109601}{40320}\leqno{(\%o2) }$$

Значение $$e$$ в форме с плавающей точкой находим, используя функцию $$float$$:

(%i3)	float(%);

$$2.71827876984127\leqno{(\%o3) }$$

Интересно провести вычисления и сравнить результаты, получающиеся для числа $$e$$ при различных степенях используемого многочлена Тейлора. Получаются следующие результаты:

$$k=1, e_{1}=1, k=2, e_{2}=2, k=3, e_{3}=2.5, k=4, e_{4}=\\ 2.666666667,k=5, e_{5}=2.708333333, k=6, e_{6}=2.716666667, e_{7}=\\ 2.718055556,k=8, e_{8}=2.718253968, k=9, e_{9}=2.718281526, e_{10}=\\ 2.718281801.$$

Отсюда видно, что значение e с точностью 0.001 вычисляется при использовании многочлена Тейлора степени не ниже 7-ой. Также следует, что число e c точностью 0.000001 или что то же самое $$10^{-6}$$ вычисляется помощи с многочлена Тейлора 9-ой или более высокой степени.

Оценку остатка ряда произведём по формуле остаточного члена ряда Маклорена $$\vert f(x_{1})-S_{n}(x_{1})\vert =\vert r_{n}(x_{1})\vert = \vert \frac{f^{n+1}(c)}{(n+1)!}\vert ,$$ где $$c$$ находится между 0 и $$x_1$$. Следует $$r_{n}(1)=\frac{e^{c}}{(n+1)!}, 0<c<1$$.. Так как $$e^{c}<e<3$$, то $$r_{n}(1)<\frac{3}{(n+1)!}$$.. При $$n = 7$$ имеем $$r_{7}<\frac{3}{7!}<0.001, e\approx 2.718$$.

Наряду с командой $$taylor$$ для разложения функций и выражений в ряды используется команда $$powerseries$$(выражение,$$x,a$$) (строится разложение для заданного выражения по переменной $$x$$ в окрестности $$a$$). Результатом выполнения команды $$powerseries$$ может быть построение её ряда Тейлора в общей форме, например:

(%i1)	powerseries(sin(x),x,0);

$$\sum_{i2=0}^{\infty }\frac{{\left( -1\right) }^{i2}\,{x}^{2\,i2+1}}{\left( 2\,i2+1\right) !}\leqno{(\%o1) }$$

(%i2)	powerseries(sin(x^2),x,0);

$$\sum_{i3=0}^{\infty }\frac{{\left( -1\right) }^{i3}\,{x}^{2\,\left( 2\,i3+1\right) }}{\left( 2\,i3+1\right) !}\leqno{(\%o2) }$$

Для разложения в ряд Тейлора функции нескольких переменных используется функция $$taylor$$ с указанием списка переменных в форме: $$taylor(expr, [x_1,x_2,... ], [a_1,a_2,... ], [n_1,n_2,... ])$$

Пример: Найти многочлен Тейлора 6-ой степени от функции $$\frac{x}{1+x}$$.

(%i1)	f(x):=x/(1+x);

$$f\left( x\right) :=\frac{x}{1+x}\leqno{(\%o1) }$$

(%i2)	powerseries(f(x),x,0);

$$x\,\sum_{i1=0}^{\infty }{\left( -1\right) }^{i1}\,{x}^{i1}\leqno{(\%o2) }$$

(%i3)	taylor(f(x),x,0,6);

$$x-{x}^{2}+{x}^{3}-{x}^{4}+{x}^{5}-{x}^{6}+\dots\leqno{(\%o3) }$$

Пример: Найти разложение функции $$\arccos(x)$$ в ряд Маклорена.

(%i6)	taylor(acos(x),x,0,12);

$$\frac{\pi }{2}-x-\frac{{x}^{3}}{6}-\frac{3\,{x}^{5}}{40}-\frac{5\,{x}^{7}}{112}-\frac{35\,{x}^{9}}{1152}-\frac{63\,{x}^{11}}{2816}+\dots\leqno{(\%o6) }$$

Пример: Найти разложение функции $$\exp(x)+1$$ по формуле Тейлора 5-ой степени в окрестности точки $$x = 2$$.

(%i7)	taylor(exp(x)+1,x,2,5);

$$1+{e}^{2}+{e}^{2}\,\left( x-2\right) +\frac{{e}^{2}\,{\left( x-2\right) }^{2}}{2}+\frac{{e}^{2}\,{\left( x-2\right) }^{3}}{6}+\frac{{e}^{2}\,{\left( x-2\right) }^{4}}{24}+\frac{{e}^{2}\,{\left( x-2\right) }^{5}}{120}+\dots\leqno{(\%o7) }$$

Пример: Найти разложение гиперболического косинуса в ряд Маклорена 8-ой степени.

taylor(cosh(x),x,10);

Получаем

$$1+\frac{1}{2}x^{2}+\frac{1}{24}x^{4}+\frac{1}{720}x^{6}+ \frac{1}{40320}x^{8}+O(x^{10}).$$

Заметим, что у аналитических функций их разложения в ряд Тейлора существуют всегда. Приведём пример функции, не имеющей разложения в ряд Тейлора и для которой команда $$taylor$$ не даёт результата: $$f(x) = 1/x^2 + x$$.

(%i8)	taylor(1/x^{2}+x,x,0,7);

$$\frac{1}{{x}^{{2}}}+x+\dots\leqno{(\%o8) }$$

В результате выполнения команды $$taylor$$ или $$powerseries$$ получаем исходное выражение $$x^{-2} + x$$. В то же время в окрестности других точек, например точки $$x = 2$$, формула Тейлора вычисляется

(%i13)	taylor(1/x^{2}+x,x,2,2);

$$\frac{2\,{2}^{{2}}+1}{{2}^{{2}}}-\frac{\left( {2}-2\,{2}^{{2}}\right) \,\left( x-2\right) }{2\,{2}^{{2}}}+\frac{\left( {{2}}^{2}+{2}\right) \,{\left( x-2\right) }^{2}}{8\,{2}^{{2}}}+\dots\leqno{(\%o13) }$$

(%i14)	ratsimp(%);

$${2}^{-{2}-3}\,\left( \left( {{2}}^{2}+{2}\right) \,{x}^{2}+\left( {2}^{{2}+3}-4\,{{2}}^{2}-8\,{2}\right) \,x+4\,{{2}}^{2}+12\,{2}+8\right) \leqno{(\%o14) }$$

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

Рассмотрим примеры такой графической визуализации для функции $$\cos(x)$$. Сравним графики самой функции $$\cos(x)$$ с графиками её разложений Тейлора различных степеней.

Пример: Сравним функцию $$\cos(x)$$ c её разложением Маклорена 4-ой степени на интервале [-5, 5].

Построим разложение

(%i15)	appr:taylor(cos(x),x,0,5);

$$1-\frac{{x}^{2}}{2}+\frac{{x}^{4}}{24}+\dots\leqno{(\%o15) }$$ (рис 3.14) Сопоставление разложения в ряд Маклорена и функции y=cos(x)

Построим график (экранная форма, в формате wxMaxima)

(%i16)	wxplot2d([appr,cos(x)], [x,-5,5], [y,-1.1,1.1],
[nticks,100]);

Выведем график в файл:

(%i17)	plot2d([appr,cos(x)], [x,-5,5], [y,-1.1,1.1],
[gnuplot_preamble, "set grid;"], [gnuplot_term, ps],
[gnuplot_out_file, "appr.eps"])$

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

Пример: Сравним функцию $$\cos(x)$$ с её разложением Маклорена 8-ой степени на интервале [-5, 5]. Сопоставим результат с предыдущим примером.

Построим разложение более высокой степени:

(%i18)	appr1:taylor(cos(x),x,0, 9);

(рис 3.15) Сопоставление двух разложений в ряд Маклорена и функции y=cos(x) $$1-\frac{{x}^{2}}{2}+\frac{{x}^{4}}{24}-\frac{{x}^{6}}{720}+\frac{{x}^{8}}{40320}+\dots\leqno{(\%o18) }$$

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

Разложение в ряд Тейлора может использоваться и для вычисления пределов (функция $$tlimit$$, по синтаксису аналогичная $$limit$$).

3.6.2 Приближённое вычисление определённых интегралов

Степенные ряды эффективны и удобны при приближённом вычислении определённых интегралов, не выражающихся в конечном виде через элементарные функции. Для вычисления $$\int_{0}^{x}f(t)dt$$ подинтегральная функция $$f(t)$$ раскладывается в степенной ряд. Если $$f(x)=a_{0}+a_{1}x+a_{2}x^{2}+\dots+a_{n}x^{n}+\dots, \ \vert x\vert<R,$$ то при $$|x| < R$$ степенной ряд можно интегрировать почленно. Получаем метод вычисления интеграла $$\int_{0}^{x}f(t)dt$$ с любой наперёд заданной точностью $$\int_{0}^{x}f(t)dt=a_{0}x+a_{1}\frac{x^{2}}{2}+ a_{2}\frac{x^{3}}{3}+\dots+a_{n}\frac{x^{n+1}}{n+1}+\dots.$$

Пример: Приближённое вычисление интеграла вероятностей

$$\phi(x)=\frac{1}{\sqrt{2\pi}}\int_{-x}^{x}e^{-t^{2}/2}dt= \frac{2}{\sqrt 2\pi }\int_{0}^{x}e^{-t^{2}/2}dt.$$

Так как $$e^{x}=1+x+\frac{x^{2}}{2!}+\frac{x^{3}}{3!}+\dots, \ \vert x\vert <\infty ,$$ то $$e^{-x^{2}/2}=1-\frac{x^{2}}{2}+ \frac{x^{4}}{2^{2}2!}-\frac{x^{6}}{2^{3}3!}+\dots$$

Подставив этот ряд под знак интеграла и произведя почленное интегрирование получаем $$\phi(x)=\frac{2}{\sqrt{2\pi}}\int_{0}^{x}e^{-t^{2}/2}dt= \frac{2}{\sqrt{2\pi}}\left[x-\frac{x^{3}}{3\cdot 2}+\frac{x^{5}} {5\cdot 2^{2}\cdot 2!}-\frac{x^{7}}{7\cdot 2^{3}\cdot 3!}+\dots\right]$$

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

Рассмотрим пример приближённого представления интеграла в виде полинома некоторой степени в том случае, когда он не вычисляется в замкнутой аналитической форме.

Пример: Вычислить $$\displaystyle \int_{0}^{1}e^{-x^{2}/2}dx$$, оценить достигнутую точность

Используем разложение подинтегральной функции в ряд. Подставляя в полученное выражение $$x = 1$$, вычисляем искомый интеграл. Так как исследуемый ряд знакопеременный, погрешность замены бесконечной суммы конечным выражением по абсолютной величине не превышает первого отброшенного члена.

(%i1)	f(x):=exp(-x^2/2);

$$f\left( x\right) :=exp\left( \frac{-{x}^{2}}{2}\right) \leqno{(\%o1) }$$

(%i2)	taylor(f(x),x,0,8);

$$1-\frac{{x}^{2}}{2}+\frac{{x}^{4}}{8}-\frac{{x}^{6}}{48}+\frac{{x}^{8}}{384}+\dots\leqno{(\%o2) }$$

Интегрируя в пределах от 0 до 1, получаем числовой результат:

(%i3) 	integrate(%,x,0,1);

$$\frac{103499}{120960}\leqno{(\%o3) }$$

(%i4)	float(%);

$$0.85564649470899\leqno{(\%o4) }$$

Точность расчёта оцениваем, интегрируя в пределах от 0 до $$a$$:

(%i5)	integrate(%o2,x,0,a);

$$\frac{35\,{a}^{9}-360\,{a}^{7}+3024\,{a}^{5}-20160\,{a}^{3}+120960\,a}{120960}\leqno{(\%o5) }$$

При $$a = 1$$ находим:

(%i6)	expand(%);

$$\frac{{a}^{9}}{3456}-\frac{{a}^{7}}{336}+\frac{{a}^{5}}{40}-\frac{{a}^{3}}{6}+a\leqno{(\%o6) }$$

(%i7)	float(1/3456);

$$2.8935185185185184\,{10}^{-4}\leqno{(\%o7) }$$

Таким образом, точность расчёта значения интеграла $$\displaystyle\int_{0}^{1}e^{-x^{2}/2}dx$$, не хуже 0,0003. Окончательно

$$\displaystyle\int_{0}^{1}e^{-x^{2}/2}dx = 2.8935 \pm 0.0003$$

3.7 Преобразование степенных рядов

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

Аппроксимацией Паде для функции $$f(x) = \sum_{k=0}^{\infty}a_k x^k$$, заданной степенным рядом, называется такая дробно-рациональная функция $$R(x) =\frac{\sum_{k=0}^{L}p_k x^k}{1+\sum_{k=1}^{M}q_k x^k}$$, чьё разложение в степенной ряд совпадает со степенным рядом $$f(x)$$ с точностью до коэффициента при $$x^{L+M}$$.

Паде-аппроксиммант задаётся значением функции в заданной точке и $$M+L$$ значениями её производных в этой же точке. Эта же информация может послужить основой для степенного ряда, так в чём же отличие? Главное отличие в том, что задав $$M+L+1$$ член степенного ряда, мы отбрасываем остальные члены ряда, приравнивая их к нулю. Паде-аппроксимант не является полиномом, поэтому задав $$M+L+1$$ членов разложения Паде-аппроксиманта в степенной ряд, мы в неявной форме задаём и остальные члены.

Чему эти дополнительные члены будут равны? Это вопрос, на который нет однозначного ответа. В одних случаях они позволят нам построить более точную аппроксимацию, в других — наоборот могут ухудшить положение. Нет способа, который позволил бы сказать, насколько точна окажется Паде-аппроксимация и в какой окрестности и с какой точностью можно получить результаты.

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

Паде-аппрокимация наиболее эффективна для функций, имеющих полюса на комплексной плоскости в окрестностях точки разложения. Например, функция $$f(x) = \frac{1}{1+\sin^2(x)}$$ непрерывна на действительной оси, но имеет полюса на комплексной плоскости. Поэтому она неэффективно аппроксимируется степенным рядом (до шестой степени включительно), но хорошо аппроксимируется по Паде со степенями числителя и знаменателя равными 4 и 2.

Функция $$pade$$, представленная в пакете Maxima, аппроксимирует отрезок ряда Тейлора, содержащий слагаемые до $$n$$-го порядка включительно, дробно-рациональной функцией. Её аргументы — ряд Тейлора, порядок числителя n, порядок знаменателя $$m$$. Разумеется, количество известных коэффициентов ряда Тейлора должно совпадать с общим количеством коэффициентов в дробно-рациональной функции минус один, поскольку числитель и знаменатель определены с точностью до общего множителя.

Синтаксис вызова функции $$pade$$:

pade (ряд Тейлора, степень числителя, степень знаменателя)

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

Пример:

(%i1)	t:taylor(exp(x),x,0,3);

$$1+x+\frac{{x}^{2}}{2}+\frac{{x}^{3}}{6}+\dots\leqno{(\%o1) }$$

(%i2)	pade(t,1,2);

$$[\frac{2\,x+6}{{x}^{2}-4\,x+6}]\leqno{(\%o2) }$$

(%i3)	taylor(sin(x)/x,x,0,7);

$$1-\frac{{x}^{2}}{6}+\frac{{x}^{4}}{120}-\frac{{x}^{6}}{5040}+\dots\leqno{(\%o3) }$$

(%i4)	pade(%,2,4);

$$[-\frac{620\,{x}^{2}-5880}{11\,{x}^{4}+360\,{x}^{2}+5880}]\leqno{(\%o4) }$$

(%i5)	taylor (1/(cos(x) - sec(x))^3, x, 0, 5);

$$-\frac{1}{{x}^{6}}+\frac{1}{2\,{x}^{4}}+\frac{11}{120\,{x}^{2}}-\frac{347}{15120}-\frac{6767\,{x}^{2}}{604800}-\frac{15377\,{x}^{4}}{7983360}+\dots\leqno{(\%o5) }$$

(%i6)	pade(%,3,inf);

$$[-\frac{120}{41\,{x}^{10}+60\,{x}^{8}+120\,{x}^{6}},\frac{8806092\,{x}^{2}-16847160}{1353067\,{x}^{10}-382512\,{x}^{8}+16847160\,{x}^{6}}]\leqno{(\%o6) }$$

Более специфичной является функция $$cf$$, которая рассчитывает коэффициенты цепной дроби, аппроксимирующей заданное выражение. Синтаксис вызова — $$cf(expr)$$. Выражение $$expr$$ должно состоять из целых чисел, квадратных корней целых чисел и знаков арифметических операций. Функция возвращает список коэффициентов (непрерывная дробь $$a + 1/(b + 1/(c + ...))$$ представляется списком [$$a,b,c,...$$]). Флаг $$cflength$$ определяет количество периодов цепной дроби. Изначально установлено значение 1. Функция $$cfdisrep$$ преобразует список (как правило выдачу функции $$''cf$$) в собственно цепную дробь вида $$a + 1/(b + 1/(c + ...)$$).

Примеры использования функций $$cf$$ и $$cfdisrep$$:

(%i1)	cf ([1, 2, -3] + [1, -2, 1]);

$$[1,1,1,2]\leqno{(\%o1) }$$

(%i2)	cfdisrep (%);

$$1+\frac{1}{1+\frac{1}{1+\frac{1}{2}}}\leqno{(\%o2) }$$

(%i3)	cflength: 3;

$$3\leqno{(\%o3) }$$

(%i4)	cf (sqrt (3));

$$[1,1,2,1,2,1,2]\leqno{(\%o4) }$$

(%i5)	cfdisrep(%);

$$1+\frac{1}{1+\frac{1}{2+\frac{1}{1+\frac{1}{2+\frac{1}{1+\frac{1}{2}}}}}}\leqno{(\%o5) }$$

(%i6)	ev (%, numer);

$$1.731707317073171\leqno{(\%o6) }$$

3.8 Решение дифференциальных уравнений в Maxima

3.8.1 Основные определения

Дифференциальным уравнением называется уравнение вида $$F(x,y,y',... ,y^{(n)}) = 0$$, где $$F(t_0,t_1,... ,t_{n+1})$$ — функция, определенная в некоторой области $$D$$ пространства $$R^{n+2}, x$$ — независимая переменная, $$y$$ — функция от $$x, y',... ,y^{(n)}$$ — ее производные.

Порядком уравнения n называется наивысший из порядков производных $$y$$, входящих в уравнение. Функция $$f(x)$$ называется решением дифференциального уравнения на промежутке $$(a; b)$$, если для всех $$x$$ из $$(a; b)$$ выполняется равенство: $$F(x,f(x),f'(x),... ,f^n(x)) = 0$$.

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

Если дифференциальное уравнение имеет одну независимую переменную, то оно называется обыкновенным дифференциальным уравнением, если же независимых переменных две или более, то такое дифференциальное уравнение называется дифференциальным уравнением в частных производных.

Пример: решить уравнение $$y' = 0$$.

Очевидно, что его решение $$f(x) = const$$ определено на $$(-\infty, \infty)$$. Отметим, что эта постоянная –– произвольная и решение — не единственное, а имеется бесконечное множество решений.

Пример: Решить уравнение $$y'=\frac{y}{x}$$, или $$\frac{dy}{dx} = \frac{y}{x}$$.

Преобразуя уравнение, получим: $$\frac{dy}{y} = \frac{dx}{x}$$. Интегрируя обе части уравнения, получим: $$\int\frac{dy}{y} = \int\frac{dx}{x} \Rightarrow ln\,y = ln \, x + ln \, C$$, или $$y = Cx$$. Общее решение изображается серией линейных интегральных кривых, проходящих через точку (0,0). При этом через любую точку, не принадлежащую (0, 0), проходит только одна интегральная кривая (решение).

Общее решение — множество решений дифференциального уравнения $$y' = f(x,y)$$ есть совокупность функций $$F(x,y,C) = 0,\, C =const$$. Частное решение получают при подстановке конкретного значения константы в общее решение. Особые решения не входят в общие решения, и через каждую точку особого решения проходит более одной интегральной кривой. Особые решения нельзя получить из общего решения ни при каких значениях постоянной $$С$$. Если построить семейство интегральных кривых дифференциального уравнения, то особое решение будет изображаться линией, которая в каждой своей точке касается по крайней мере одной интегральной кривой.

Пример: Рассмотрим уравнение $$y'=\frac{-x}{y}$$. Преобразуя его, найдём: $$\frac{dy}{dx} = - \frac{x}{y} \Rightarrow 2ydy + 2xdx = 0 \Leftrightarrow d(x^2+y^2) = 0$$. Интегрируя, получаем $$x^2 + y^2 = C$$.

Пример: Дифференциальное уравнение $$y' = 2 \sqrt{y}$$ имеет общее решение $$y = (x-C)^2$$ и особое решение $$y = 0$$. При конкретном значении $$С$$ (например, $$C = 1$$) получаем частное решение: $$y = (x - 1)^2$$.

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

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

Пусть $$f(x,y)$$ — непрерывная функция в области $$D = {(x,y; a < x < b; c < y < d)}$$, причём частная производная $$\frac{\partial f}{\partial y}(x,y)$$ также непрерывна в $$D$$. Тогда существует единственное решение $$y = y(x)$$ дифференциального уравнения $$y' = f(x,y)$$ с начальным условием $$y(x_0) = y_0,(x_0,y_0) \in D$$. Следовательно, через точку $$(x_0,y_0) \in D$$ проходит только одна интегральная кривая.

3.8.2 Функции для решения дифференциальных уравнений в Maxima

В Maxima предусмотрены специальные средства решения задачи Коши для систем обыкновенных дифференциальных уравнений, заданных как в явной форме $$\frac{dx}{dt}= F(t,x)$$, так и в неявной $$M\frac{dy}{dt}=F(t,x)$$, где $$M$$ — матрица, — т.н. решатель ОДУ ($$solver ODE$$), обеспечивающий пользователю возможность выбора метода, задания начальных условий и др.

Функция $$ode2$$ позволяет решить обыкновенные дифференциальные уравнения первого и второго порядков.

Синтаксис вызова $$ode2(eqn,dvar,ivar)$$, где $$eqn$$ — выражение, определяющее само дифференциальное уравнение, зависимая переменная — $$dvar$$, независимая переменная — $$ivar$$.

Дифференциальное уравнение представляется в форме с "замороженной" производной (т.е. с производной, вычисление которой запрещено с помощью одиночной кавычки: $$'diff(y,x)$$). Другой вариант явно указать зависимость $$y = y(x)$$ — использовать функцию $$depends$$ (в этом случае можно не использовать начальный апостроф см. при- мер). Если ode2 не может получить решение, она возвращает значение $$false$$.

Посредством функции $$ode2$$ могут быть решены следующие типы ОДУ первого порядка: линейные, ОДУ с разделяющимися переменными, однородные ОДУ, уравнения в полных дифференциалах, уравнения Бернулли, обобщённые однородные уравнения.

Кроме того, при помощи функции $$ode2$$ могут быть решены следующие типы уравнений второго порядка: с постоянными коэффициентами; в полных дифференциалах; линейные однородные с переменными коэффициентами, которые могут быть сведены к уравнениям с постоянными коэффициентами; уравнения Эйлера; уравнения, разрешимые методом вариации постоянных; уравнения, свободные от независимой переменной, допускающие понижение порядка.

Тип используемого метода сохраняется в переменной $$method$$. При использовании интегрирующего множителя он сохраняется в переменной $$intfactor$$. Частное решение неоднородного уравнения сохраняется в переменной $$yp$$.

Для отыскания частных решений задач Коши с начальными условиями используются функции $$ic1$$ (для уравнений первого порядка) и $$ic2$$ (для уравнений второго порядка). Частные решения граничных задач для уравнений второго порядка используют функцию $$bc2$$.

Рассмотрим примеры использования функции $$ode2$$.

Вариант использования отложенного дифференцирования ($$'diff$$):

(%i1)	ode2('diff(y,x)=2*y+exp(x),y,x);

$$y=\left( \%c-{e}^{-x}\right) \,{e}^{2\,x}\leqno{(\%o1) }$$

Вариант с явным указанием зависимости $$y = y(x)$$:

(%i1)	depends(y,x);

$$[y\left( x\right) ]\leqno{(\%o1) }$$

(%i2)	ode2(diff(y,x)=2*y+exp(x),y,x);

$$y=\left( \%c-{e}^{-x}\right) \,{e}^{2\,x}\leqno{(\%o2) }$$

Параметр $$\%с$$ — постоянная интегрирования для уравнения первого порядка.

Решение уравнения второго порядка:

(%i4)	ode2('diff(y,x,2)-3*'diff(y,x)+2*y=0,y,x);

$$y=\%k1\,{e}^{2\,x}+\%k2\,{e}^{x}\leqno{(\%o4) }$$

Параметры $$\%k1$$ и $$\%k2$$ — постоянные интегрирования для уравнений второго порядка.

Рассмотрим варианты вычисления частных решений: для уравнения первого порядка

(%i5)	ic1(%o1,x=1,y=1);

$$y={e}^{-2}\,\left( \left( e+1\right) \,{e}^{2\,x}-{e}^{x+2}\right) \leqno{(\%o5) }$$

для уравнения второго порядка

(%i6)	ic2(%o4,x=0,y=1,diff(y,x)=1);

$$y={e}^{x}\leqno{(\%o6) }$$

3.8.3 Решение основных типов дифференциальных уравнений

3.8.3.1 Уравнения с разделяющимися переменными

Уравнениями с разделяющимися переменными называются уравнения вида $$y' = f(x)\cdot g(y)$$, где $$f(x)$$ — функция, непрерывна на некотором интервале ($$a,b$$), а функция $$g(y)$$ — функция, непрерывна на интервале ($$c,d$$), причем $$g(y) \ne 0$$ на ($$c,d$$)).

Преобразуя уравнение, получаем: $$\frac{dy}{dx}=f(x)\cdot g(y)\Leftrightarrow \frac{dy}{g(y)}=f(x)dx$$

Интегрируя обе части, получаем $$\int\frac{dy}{g(y)}=\int f(x)dx$$. Обозначая $$G(y)$$ любую первообразную для $$\frac{1}{g(y)}$$ , а $$F(x)$$ — любую первообразную для $$f(x)$$, получаем общее решение дифференциального уравнения в виде неявно выраженной функции $$G(y) = F(x) + C$$.

Пример решения в Maxima:

Отыскиваем общее решение:

(%i1) difur1:'diff(y,x)=sqrt(1-y^2)/sqrt(1-x^2);

$$\frac{d}{d\,x}\,y=\frac{\sqrt{1-{y}^{2}}}{\sqrt{1-{x}^{2}}}\leqno{(\%o1) }$$

(%i2) rez:ode2(difur1,y,x);

$$asin\left( y\right) =asin\left( x\right) +\%c\leqno{(\%o2) }$$

Отыскиваем различные варианты частных решений:

(%i3) ic1(rez,x=0,y=0);

$$asin\left( y\right) =asin\left( x\right) \leqno{(\%o3) }$$

(%i4) ic1(rez,x=0,y=1);

$$asin\left( y\right) =\frac{2\,asin\left( x\right) +\pi }{2}\leqno{(\%o4) }$$

3.8.3.2 Однородные уравнения

Под однородными уравнениями понимаются уравнения вида $$y'= f(\frac{y}{x})$$. Для их решения используется замена вида $$y = u \cdot x$$, после подстановки которой получается уравнение с разделяющимися переменными: $$y' = u'\,x+u \! \Rightarrow \! u'\,x+u = f(u)$$. Разделяя переменные и интегрируя, получаем: $$x\frac{du}{dx}=f(u)-u \Rightarrow \int\frac{du}{f(u)-u}\,=\, \int\frac{dx}{x}$$

Пример решения в Maxima:

Находим общее решение:

(%i1)	homode:'diff(y,x) = (y/x)^2 + 2*(y/x);

$$\frac{d}{d\,x}\,y=\frac{{y}^{2}}{{x}^{2}}+\frac{2\,y}{x}\leqno{(\%o1) }$$

(%i2)	ode2(homode,y,x);

$$-\frac{x\,y+{x}^{2}}{y}=\%c\leqno{(\%o2) }$$

Находим частное решение:

(%i3)	ic1(%,x=2,y=1);

$$-\frac{x\,y+{x}^{2}}{y}=-6\leqno{(\%o3) }$$

Более общий вариант дифференциальных уравнений, уравнения вида: $$y'= \frac{a_1x+b_1y+c_1}{a_2x+b_2y+c_2}$$ — сводим их к однородным. Maxima не способна решать такие уравнения при помощи $$ode2$$ непосредственно, а лишь после необходимого преобразования.

3.8.3.3 Линейные уравнения первого порядка

Дифференциальное уравнение называется линейным относительно неизвестной функции и ее производной, если оно может быть записано в виде: $$y' = P(x)\cdot y = Q(x)$$ при этом, если правая часть $$Q(x)$$ равна нулю, то такое уравнение называется линейным однородным дифференциальным уравнением, если правая часть $$Q(x)$$ не равна нулю, то такое уравнение называется линейным неоднородным дифференциальным уравнением. При этом $$P(x)$$ и $$Q(x)$$ — функции непрерывные на некотором промежутке $$x \in (a,b)$$.

Рассмотрим решение линейного дифференциального уравнения в Maxima:

(%i1)	lineq1:'diff(y,x)-y/x=x;

$$\frac{d}{d\,x}\,y-\frac{y}{x}=x\leqno{(\%o1) }$$

(%i2)	ode2(lineq1,y,x);

$$y=x\,\left( x+\%c\right) \leqno{(\%o2) }$$

При работе с Maxima не требуется приводить дифференциальное уравнение к стандартной форме вида

$$y' = P(x)\cdot y = Q(x)$$

(%i3)	lineq2:y^2-(2*x*y+3)*'diff(y,x)=0;

$${y}^{2}-\left( 2\,x\,y+3\right) \,\left( \frac{d}{d\,x}\,y\right) =0\leqno{(\%o3) }$$

(%i4)	ode2(lineq2,y,x);

$$\frac{x\,y+1}{{y}^{3}}=\%c\leqno{(\%o4) }$$

3.8.3.4 Уравнения в полных дифференциалах

Дифференциальное уравнение первого порядка вида: $$P(x,y)dx+Q(x,y)dy = 0$$ называется уравнением в полных дифференциалах, если левая часть этого уравнения представляет собой полный дифференциал некоторой функции $$u = F(x,y)$$. Данное дифференциальное уравнение является уравнением в полных дифференциалах, если выполняется условие:

$$\frac{\partial P}{\partial y} = \frac{\partial Q}{\partial x}.$$

Общий интеграл уравнения имеет вид $$U(x,y) = 0$$.

Если уравнение $$P(x,y)dx + Q(x,y)dy =$$0 не является уравнением в полных дифференциалах, но выполняются условия теоремы единственности, то существует функция $$\mu = \mu (x,y)$$ (интегрирующий множитель) такая, что

$$\mu (Pdx + Qdy) = dU$$

Функция $$\mu$$ удовлетворяет условию:

$$\frac{\partial (\mu P)}{\partial y} = \frac{\partial (\mu Q)}{\partial x}$$

Примеры решения в Maxima:

Для решения в Maxima дифференциальное уравнение представляется в форме

$$P(x,y)+Q(x,y) \frac{dy}{dx} = 0$$

Уравнение, приводимое к уравнению в полных дифференциалах

(%i1)	deq:(2*x*y+x^2*y+y^3/3)+(x^2+y^2)*'diff(y,x)=0;

$$\left( {y}^{2}+{x}^{2}\right) \,\left( \frac{d}{d\,x}\,y\right) +\frac{{y}^{3}}{3}+{x}^{2}\,y+2\,x\,y=0\leqno{(\%o1) }$$

(%i2)	ode2(deq,y,x);

$$\frac{{e}^{x}\,{y}^{3}+3\,{x}^{2}\,{e}^{x}\,y}{3}=\%c\leqno{(\%o2) }$$

Указание на интегрирующий множитель

(%i3)	intfactor;

$${e}^{x}\leqno{(\%o3) }$$

Указание на использованный метод

(%i4)	method;

$$exact\leqno{(\%o4) }$$

Уравнение в полных дифференциалах

(%i5)	deq1:(3*x^2+6*x*y^2)+(6*x^2*y+4*y^3)*'diff(y,x)=0;

$$\left( 4\,{y}^{3}+6\,{x}^{2}\,y\right) \,\left( \frac{d}{d\,x}\,y\right) +6\,x\,{y}^{2}+3\,{x}^{2}=0\leqno{(\%o5) }$$

(%i6)	ode2(deq1,y,x);

$${y}^{4}+3\,{x}^{2}\,{y}^{2}+{x}^{3}=\%c\leqno{(\%o6) }$$

Указание на использованный метод

(%i7)	method;

$$exact\leqno{(\%o7) }$$

3.8.3.5 Уравнения Бернулли

Уравнением Бернулли называется уравнение вида $$y' + P(x)\cdot y = Q(x) \cdot y^{\alpha}$$ где $$P$$ и $$Q$$ — функции от x или постоянные числа, а $$\alpha$$ — постоянное число, не равное 0 и 1При $$\alpha = 0$$ получаем неоднородное, а при $$\alpha = 1$$ однородное линейное уравнение..

Для решения уравнения Бернулли применяют подстановку $$z=\frac{1}{y^{\alpha-1}}$$ , с помощью которой, уравнение Бернулли приводится к линейному.

Пример решения уравнения Бернулли с помощью Maxima:

(%i1)	deq:'diff(y,x)=4/x*y+x*sqrt(y);

$$\frac{d}{d\,x}\,y=\frac{4\,y}{x}+x\,\sqrt{y}\leqno{(\%o1) }$$

(%i2)	ode2(deq,y,x);

$$y={x}^{4}\,{\left( \frac{log\left( x\right) }{2}+\%c\right) }^{2}\leqno{(\%o2) }$$

(%i3)	method;

$$bernoulli\leqno{(\%o3) }$$

(%i4)	de1:'diff(y,x)+y/x=-x*y^2;

$$\frac{d}{d\,x}\,y+\frac{y}{x}=-x\,{y}^{2}\leqno{(\%o4) }$$

(%i5)	ode2(de1,y,x);

$$y=\frac{1}{x\,\left( x+\%c\right) }\leqno{(\%o5) }$$

3.8.3.6 Уравнения высших порядков

В Maxima при помощи функции $$ode2$$ возможно прямое решение лишь линейных дифференциальных уравнений второго порядка При решении выполняется проверка, является ли заданное уравнение линейным, т.е. возможно ли его преобразование к форме $$y'' + p(x)y + q(x)y = r(x)$$.

Первоначально отыскивается решение однородного уравнения вида $$y'' + p(x)y + q(x)y = 0$$ в форме $$y = k_1 y_1 + k_2 y_2$$ ($$k_1,k_2$$ — произвольные постоянные). Если $$r(x) \ne 0$$, отыскивается частное решение неоднородного уравнения методом вариации постоянных. . .

3.8.3.7 Уравнения с постоянными коэффициентами

Решения однородных уравнений вида $$y'' + a * y' + b * y = 0$$ отыскиваются по результатам решения характеристического уравнения $$r^2 + ar + b = 0$$. Возможны следующие варианты комбинаций его корней $$r_1, r_2$$:

$$r_1, r_2$$ — вещественные и различные. Решение представляется в форме $$y=k_1 \cdot e^{r_1 \cdot x} + k_2 \cdot e^{r_2 \cdot x}$$.

$$r_1, = r_2$$ — корни вещественные одинаковые. Решение представляется в форме $$y= (k_1 + k_2 \cdot x) e^{r_1 \cdot x}$$.

$$r_1, r_2$$ — комплексные (сопряжённые). Если $$r_1 = \alpha + \beta i ,\ r_2 = \alpha - \beta i$$, то решение представляется в виде $$y = e^{\alpha x}(k_1 \cos \, \beta x + k_2 \sin \beta x)$$

Общее решение неоднородного уравнения с постоянными коэффициентами представляется в виде суммы общего решения соответствующего однородного уравнения и какого-либо частного решения неоднородного.

Примеры решения ОДУ второго порядка с постоянными коэффициентами:

Неоднородное уравнение общего вида:

(%i1)	de1:2*'diff(y,x,2)-'diff(y,x)-y=4*x*exp(2*x);

$$2\,\left( \frac{{d}^{2}}{d\,{x}^{2}}\,y\right) -\frac{d}{d\,x}\,y-y=4\,x\,{e}^{2\,x}\leqno{(\%o1) }$$

(%i2)	ode2(de1,y,x);

$$y=\frac{\left( 20\,x-28\right) \,{e}^{2\,x}}{25}+\%k1\,{e}^{x}+\%k2\,{e}^{-\frac{x}{2}}\leqno{(\%o2) }$$

Частное решение неоднородного уравнения сохраняется в переменной $$yp$$:

(%i3)	yp;

$$\frac{\left( 20\,x-28\right) \,{e}^{2\,x}}{25}\leqno{(\%o3) }$$

Неоднородное уравнение с кратными корнями характеристического уравнения:

(%i1)	de2:'diff(y,x,2)-2*'diff(y,x)+y=x*exp(x);

$$\frac{{d}^{2}}{d\,{x}^{2}}\,y-2\,\left( \frac{d}{d\,x}\,y\right) +y=x\,{e}^{x}\leqno{(\%o1) }$$

(%i2)	ode2(de2,y,x);

$$y=\frac{{x}^{3}\,{e}^{x}}{6}+\left( \%k2\,x+\%k1\right) \,{e}^{x}\leqno{(\%o2) }$$

(%i3)	yp;

$$\frac{{x}^{3}\,{e}^{x}}{6}\leqno{(\%o3) }$$

Неоднородное уравнение с комплексными корням:

(%i4)	de3:'diff(y,x,2)+y=x*sin(x);

$$\frac{{d}^{2}}{d\,{x}^{2}}\,y+y=x\,sin\left( x\right) \leqno{(\%o4) }$$

(%i5)	ode2(de3,y,x);

$$y=\frac{2\,x\,sin\left( x\right) +\left( 1-2\,{x}^{2}\right) \,cos\left( x\right) }{8}+\%k1\,sin\left( x\right) +\%k2\,cos\left( x\right) \leqno{(\%o5) }$$

(%i6)	yp;

$$\frac{2\,x\,sin\left( x\right) +\left( 1-2\,{x}^{2}\right) \,cos\left( x\right) }{8}\leqno{(\%o6) }$$

3.8.3.8 Уравнения с переменными коэффициентами

Аналогично уравнению с постоянными коэффициентами, общее решение однородного уравнения $$y''+p(x)y'+q(x)y=0$$ имеет вид $$y=C_1 y_1 + C_2 y_2$$, где $$y_1,y_2$$ — линейно независимые решения однородного ОДУ (фундаментальная система решений).

Общее решение неоднородного уравнения $$y''+p(x)y'+q(x)y=f(x)$$ с непрерывными коэффициентами и правой частью имеет вид y$$ = y_0 + Y$$ , где $$y_0$$ — общее решение соответствующего однородного уравнения, $$Y$$ — частное решение неоднородного.

Если известна фундаментальная система решений однородного уравнения, общее решение неоднородного может быть представлено в форме: $$y = C_1(x) \cdot y_1 + C_2(x) \cdot y_2 ,$$ где $$C_1(x),C_2(x)$$ определяются методом вариации произвольных постоянных.

Пример:

(%i3)	difur:x^2*'diff(y,x,2)-x*'diff(y,x)=3*x^3;

$${x}^{2}\,\left( \frac{{d}^{2}}{d\,{x}^{2}}\,y\right) -x\,\left( \frac{d}{d\,x}\,y\right) =3\,{x}^{3}\leqno{(\%o3) }$$

(%i4)	ode2(difur,y,x);

$$y={x}^{3}+\%k2\,{x}^{2}-\frac{\%k1}{2}\leqno{(\%o4) }$$

Пример:

(%i3)	difur1:x*'diff(y,x,2)+'diff(y,x)=x^2;

$$x\,\left( \frac{{d}^{2}}{d\,{x}^{2}}\,y\right) +\frac{d}{d\,x}\,y={x}^{2}\leqno{(\%o3) }$$

(%i4)	ode2(difur1,y,x);

$$y=\%k1\,log\left( x\right) +\frac{{x}^{3}}{9}+\%k2\leqno{(\%o4) }$$

3.8.3.9 Уравнение Эйлера

Однородное уравнение $$x^2y''+axy'+by = 0$$ называется уравнением Эйлера. Его общее решение имеет вид $$y = k_1 x^{r_1} + k_2 x^{r_2}$$, где $$r_1$$ и $$r_2$$ — решения уравнения r(r _ 1) + ar + b = 0.

В случае, когда уравнение $$r(r - 1) + ar + b = 0$$ имеет двукратный корень $$r$$, решение представляется в форме $$y = k_1 x^r +k_2 ln(x) x^r$$.

Неоднородное уравнение типа Эйлера сводится к однородному с постоянными коэффициентами путём соответствующей замены.

Пример:

(%i1)	du:x^2*'diff(y,x,2)+x*'diff(y,x)+y=1;

$${x}^{2}\,\left( \frac{{d}^{2}}{d\,{x}^{2}}\,y\right) +x\,\left( \frac{d}{d\,x}\,y\right) +y=1\leqno{(\%o1) }$$

(%i2)	ode2(du,y,x);

$$y={\sin(\log(x))}^{2}+\%k1 \sin(\log(x)) +{\cos(\log(x))}^{2}+\%k2 \cos(\log(x)) \leqno{(\%o2)}$$

3.8.3.10 Граничные задачи

Для задания граничных условий при интегрировании ОДУ второго порядка используется функция $$bc2$$.

Синтаксис вызова: $$bc2 (solution, xval1, yval1, xval2, yval2)$$, где $$xval1$$ — значение $$x$$ в первой граничной точке, $$yval1$$ — значение решения $$y$$ в той же точке (обе величины задаются в форме $$x = a,y = b$$).

Пример использования $$ode2$$ и $$bc2$$:

(%i1)	'diff(y,x,2) + y*'diff(y,x)^3 = 0;

$$\frac{{d}^{2}}{d\,{x}^{2}}\,y+y\,{\left( \frac{d}{d\,x}\,y\right) }^{3}=0\leqno{(\%o1) }$$

(%i2)	ode2(%,y,x);

$$\frac{{y}^{3}+6\,\%k1\,y}{6}=x+\%k2\leqno{(\%o2) }$$

(%i3)	bc2(%,x=0,y=1,x=1,y=3);

$$\frac{{y}^{3}-10\,y}{6}=x-\frac{3}{2}\leqno{(\%o3) }$$

3.8.4 Операторный метод решения

Для решения систем обыкновенных линейных дифференциальных уравнений в Maxima имеется функция $$desolve$$. Работа функции $$desolve$$ основана на преобразовании Лапласа заданных дифференциальных уравнений.

Пусть задана функция действительного переменного $$f(t)$$, которая удовлетворяет следующим условиям:

однозначна и непрерывна вместе со своими производными $$n$$-го порядка для всех $$t > 0$$, кроме тех, где она и ее производные имеют разрывы 1-го рода. При этом в каждом конечном интервале изменения имеется конечное число точек разрыва;

$$f(t) = 0$$ для всех $$t > 0$$;

возрастает медленнее некоторой экспоненциальной функции $$M\cdot e^{at}$$, где $$M$$ и $$a$$ — некоторые положительные величины, т.е. всегда можно указать такие $$M$$ и $$a$$, чтобы при любом $$t > 0$$ соблюдалось неравенство $$|f(t)| < M \cdot e^{at}$$.

Рассматриваемой функции $$f(t)$$ ставится в соответствие новая функция, определяемая равенством $$F(s) = L \{ f (t)\} = \int \limits_0^{\infty} e^{#8722;st} f (t) dt,$$ где $$s$$ — положительное действительное число или комплексное число с положительной действительной частью.

Функция $$f(t)$$ при этом называется оригиналом, а $$F(s)$$ — изображением функции $$f(t)$$ по Лапласу. Переход от оригинала к изображению называется преобразованием Лапласа. Соответственно, обратный переход от изображения к оригиналу называется обратным преобразованием Лапласа.

Для преобразования Лапласа выполняется теорема единственности: если две непрерывные функции $$f(x)$$ и $$g(x)$$ имеют одно и то же изображение по Лапласу $$F(p)$$, то они тождественно равны.

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

При вычислении преобразования Лапласа производные заменяются алгебраическими выражениями следующего вида: $$pF(p) - f(0) = f'(t)$$ $$p^2F(p) - pf(0) - f'(0) = f''(t)$$ и т.д., поэтому использование преобразования Лапласа для решения систем ОДУ требует задания начальных условий.

Использование $$desolve$$ ограничивается одним из свойств преобразования Лапласа: если $$L{f(t)} = F(s)$$, то $$L{tf(t)} = -F(s)$$. Поэтому $$desolve$$ предполагает, что решается система ОДУ с постоянными коэффициентами.

Синтаксис вызова $$desolve: desolve(delist,fnlist)$$, где $$delist$$ — список решаемых дифференциальных уравнений, $$fnlist$$ — список искомых функций. При использовании $$desolve$$ необходимо явно задавать функциональные зависимости (вместо $$'diff(y,x)$$ использовать запись $$diff(y(x),x)$$).

Примеры использования $$desolve$$:

Система ОДУ первого порядка:

(%i1)	de1:diff(f(x),x)=diff(g(x),x)+sin(x);

$$\frac{d}{d\,x}\,f\left( x\right) =\frac{d}{d\,x}\,g\left( x\right) +sin\left( x\right) \leqno{(\%o1) }$$

(%i2)	de2:diff(g(x),x,2)=diff(f(x),x) - cos(x);

$$\frac{{d}^{2}}{d\,{x}^{2}}\,g\left( x\right) =\frac{d}{d\,x}\,f\left( x\right) -cos\left( x\right) \leqno{(\%o2) }$$

(%i3)	desolve([de1,de2],[f(x),g(x)]);

$$[f\left( x\right) ={e}^{x}\left( \left. \frac{d}{dx}\,g\left( x\right) \right|_{x=0}\right) -\left. \frac{d}{dx}g\left( x\right) \right|_{x=0}+f\left( 0\right) ,\leqno{(\%o3) }$$ $$g\left( x\right) ={e}^{x}\left( \left. \frac{d}{dx}g\left( x\right) \right|_{x=0}\right) -\left. \frac{d}{dx}g\left( x\right) \right|_{x=0}+cos\left( x\right) +g\left( 0\right) -1]$$

Единичное дифференциальное уравнение второго порядка:

(%i1)	de3:diff(f(x),x,2)+f(x) = 2*x;

$$\frac{{d}^{2}}{d\,{x}^{2}}\,f\left( x\right) +f\left( x\right) =2\,x\leqno{(\%o1) }$$

(%i2)	desolve(de3,f(x));

$$f\left( x\right) =sin\left( x\right) \,\left( \left. \frac{d}{d\,x}\,f\left( x\right) \right|_{x=0}-2\right) +f\left( 0\right) \,cos\left( x\right) +2\,x\leqno{(\%o2) }$$

Для указания начальных условий используется функция $$atvalue$$.

Синтаксис вызова: a$$tvalue(expr, [x_1 = a_1,...,x_m = a_m],c) atvalue(expr,x_1 = a_1,c)$$

Функция $$atvalue$$ присваивает значение c выражению $$expr$$ в точке $$x = a$$. Выражение $$expr$$ — функция $$f(x_1,... ,x_m)$$ или производной $$diff(f(x_1,...,x_m),x_1,n_1,...,x_n,n_m)$$ Здесь $$n_i$$ — порядок дифференцирования по переменной $$x_i$$.

Пример использования $$desolve$$ и $$atvalue$$:

(%i1)	de1:diff(f(x),x)=diff(g(x),x)+sin(x);

$$\frac{d}{d\,x}\,f\left( x\right) =\frac{d}{d\,x}\,g\left( x\right) +sin\left( x\right) \leqno{(\%o1) }$$

(%i2)	de2:diff(g(x),x,2)=diff(f(x),x) - cos(x);

$$\frac{{d}^{2}}{d\,{x}^{2}}\,g\left( x\right) =\frac{d}{d\,x}\,f\left( x\right) -cos\left( x\right) \leqno{(\%o2) }$$

(%i3)	atvalue(f(x),x=0,1);

$$1\leqno{(\%o3) }$$

(%i4)	atvalue(g(x),x=0,2);

$$2\leqno{(\%o4) }$$

(%i5)	atvalue(diff(g(x),x),x=0,3);

$$3\leqno{(\%o5) }$$

(%i6)	desolve([de1,de2],[f(x),g(x)]);

$$[f\left( x\right) =3\,{e}^{x}-2,g\left( x\right) =cos\left( x\right) +3\,{e}^{x}-2]\leqno{(\%o6) }$$

Управление начальными условиями осуществляется при помощи функций $$properties$$ и $$prontprops$$. Функция $$properties$$ (синтаксис вызова — $$properties(a)$$) печатает свойства переменной (атома $$a$$), а функция $$printprops$$ печатает информацию о заданном свойстве переменной. Кроме того, функция $$at$$ вычисляет значение выражения в заданной точке с учетом свойства $$atvalue$$.

Синтаксис вызова $$printprops$$:

$$printprops(a,i)$$

$$printprops([a1,... ,an],i)$$

$$printprops(all,i)$$

Данная функция позволяет просмотреть свойства атома $$a$$ (или группы атомов Lisp, указанных в списке), определённые индикатором $$i$$.

Отмена установок, произведённых $$atvalue$$, осуществляется функцией $$remove$$ (удаление свойства $$p$$ у атомов $$a_1,... ,a_n$$ осуществляется вызовом remove(a1,p1,... ,an,pn); удаление списка свойств — вызовом $$remove([a_1,... ,a_m], [p_1,... ,p_n],...)$$).

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

(%i1)	eq1:'diff(f(x),x)='diff(g(x),x)+sin(x);

$$\frac{d}{d\,x}\,f\left( x\right) =\frac{d}{d\,x}\,g\left( x\right) +sin\left( x\right) \leqno{(\%o1) }$$

(%i2)	eq2:'diff(g(x),x,2)='diff(f(x),x)-cos(x);

$$\frac{{d}^{2}}{d\,{x}^{2}}\,g\left( x\right) =\frac{d}{d\,x}\,f\left( x\right) -cos\left( x\right) \leqno{(\%o2) }$$

(%i3)	atvalue('diff(g(x),x),x=0,a);

$$a\leqno{(\%o3) }$$

(%i4)	atvalue(f(x),x=0,1);

$$1\leqno{(\%o4) }$$

(%i5)	properties(f);

$$[atvalue]\leqno{(\%o5) }$$

(%i6)	printprops(f,atvalue);

$$f\left( 0\right) =1$$ $$done\leqno{(\%o6) }$$

(%i7)	desolve([eq1,eq2],[f(x),g(x)]);

$$[f\left( x\right) =a\,{e}^{x}-a+1,g\left( x\right) =cos\left( x\right) +a\,{e}^{x}-a+g\left( 0\right) -1]\leqno{(\%o7) }$$

(%i8)	at(%,[x=1]);

$$[f\left( 1\right) =e\,a-a+1,g\left( 1\right) =e\,a-a+cos\left( 1\right) +g\left( 0\right) -1]\leqno{(\%o8) }$$

Ещё один пример анализа свойств:

(%i9)	atvalue (f(x,y), [x = 0, y = 1], a^2);

$${a}^{2}\leqno{(\%o9) }$$

(%i10)	atvalue ('diff (f(x,y), x), x = 0, 1 + y);

$$@2+1\leqno{(\%o10) }$$

(%i11)	printprops (all, atvalue);

$$\left. \frac{d}{d\,@1}\,g\left( @1\right) \right|_{[@1=0]}=a\\ \left. \frac{d}{d\,@1}\,f\left( @1,@2\right) \right|_{[@1=0]}=@2+1\\ f\left( 0,1\right) ={a}^{2}\\ f\left( 0\right) =1$$ $$done\leqno{(\%o11) }$$

3.8.5 Дополнительные возможности решения ОДУ

3.8.5.1 Пакет contrib_ode

Как видно из описания возможностей Maxima выше, возможности основной функции для аналитического решения ОДУ — функции ode2 — весьма ограничены. Для расширения возможностей решения ОДУ первого и второго порядка в последних версиях Maxima существует пакет расширения contrib_ode. При помощи $$contrib_ode$$ возможно решение уравнений Клеро, Лагранжа, Риккати и др. В общем случае результат — список решений. Для некоторых уравнений (в частности Риккати) решение представляется в форме другого ОДУ — результата замены переменных. Функция $$contrib_ode$$ реализует методы факторизации ($$factorization$$), Клеро (Clairault), Лагранжа (Lagrange), Риккати (Riccati), Абеля (Abel) и метод симметрии Ли (Lie symmetry method).

Для использования пакет contrib_ode необходимо загрузить:

(%i1)	load("contrib_ode")$

Пример решения ОДУ с использованием функции $$contrib_ode$$:

(%i2)	eqn:x*'diff(y,x)^2-(1+x*y)*'diff(y,x)+y=0;

$$x\,{\left( \frac{d}{d\,x}\,y\right) }^{2}-\left( x\,y+1\right) \,\left( \frac{d}{d\,x}\,y\right) +y=0\leqno{(\%o2) }$$

(%i3)	contrib_ode(eqn,y,x);

$$x\,{\left( \frac{d}{d\,x}\,y\right) }^{2}-\left( x\,y+1\right) \,\left( \frac{d}{d\,x}\,y\right) +y=0\leqno{(\%t3) }$$

first order equation not linear in y'

$$[y=log\left( x\right) +\%c,y=\%c\,{e}^{x}]\leqno{(\%o3) }$$

(%i4)	method;

$$factor\leqno{(\%o4) }$$

Достоинство $$contrib_ode$$ — возможность решения нелинейных ОДУ первого порядка, т.к. они могут иметь в общем случае несколько решений, результат представляется в виде списка.

Синтаксис вызова $$contrib_ode$$ не отличается от синтаксиса вызова $$ode2$$.

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

3.8.5.2 Уравнения Клеро и Лагранжа

Уравнение Клеро

(%i1)	load("contrib_ode")$
(%i2)	eqn:'diff(y,x)^2+x*'diff(y,x)-y=0;

$${\left( \frac{d}{d\,x}\,y\right) }^{2}+x\,\left( \frac{d}{d\,x}\,y\right) -y=0\leqno{(\%o2) }$$

(%i3)	contrib_ode(eqn,y,x);

$${\left( \frac{d}{d\,x}\,y\right) }^{2}+x\,\left( \frac{d}{d\,x}\,y\right) -y=0\leqno{(\%t3) }$$

first order equation not linear in y'

$$[y=\%c\,x+{\%c}^{2},y=-\frac{{x}^{2}}{4}]\leqno{(\%o3) }$$

(%i4)	method;

$$clairault\leqno{(\%o4) }$$

Уравнение Лагранжа

(%i5)	leq:y=(1+'diff(y,x))*x+('diff(y,x))^2;

$$y={\left( \frac{d}{d\,x}\,y\right) }^{2}+x\,\left( \frac{d}{d\,x}\,y+1\right) \leqno{(\%o5) }$$

(%i6)	contrib_ode(leq,y,x);

$$y={\left( \frac{d}{d\,x}\,y\right) }^{2}+x\,\left( \frac{d}{d\,x}\,y+1\right)\leqno{(\%t6)}$$

first order equation not linear in y'

$$[[x={e}^{-\%t}\,\left( \%c-2\,\left( \%t-1\right) \,{e}^{\%t}\right) ,y=\left( \%t+1\right) \,x+{\%t}^{2}]]\leqno{(\%o6) }$$

(%i7)	method;

$$lagrange\leqno{(\%o7) }$$

В некоторых случаях возможно только решение в параметрической форме. Пример ($$\%t$$ — параметр):

(%i8	 eqn:'diff(y,x)=(x+y)^2;

$$\frac{d}{d\,x}\,y={\left( y+x\right) }^{2}\leqno{(\%o8) }$$

(%i9)	contrib_ode(eqn,y,x);

$$[[x=\%c-atan\left( \sqrt{\%t}\right) ,y=-x-\sqrt{\%t}],\leqno{(\%o9) }$$ $$[x=atan\left( \sqrt{\%t}\right) +\%c,y=\sqrt{\%t}-x]]$$

(%i10)	method;

$$lagrange\leqno{(\%o10) }$$

3.8.5.3 Другие задачи с использованием contrib_ode

Пакет $$contrib_ode$$ позволяет решать дифференциальные уравнения, не разрешимые при помощи $$ode2$$ непосредственно. Пример — обобщённые однородные уравнения (см. выше). Представленные задачи используют методы Абеля и симметрии Ли.

(%i11)	eqn:(2*x-y+4)*'diff(y,x)+(x-2*y+5)=0;

$$\left( -y+2\,x+4\right) \,\left( \frac{d}{d\,x}\,y\right) -2\,y+x+5=0\leqno{(\%o11) }$$

(%i12)	contrib_ode(eqn,y,x);

$$\parbox{8ex}{(\%o12)}[\frac{log\left(3-\frac{2(2x+4)-x-5}{-y+2x+4}\right)-3log\left(1-\frac{2(2x+4)-x-5}{-y+2x+4}\right)+2log\left(-\frac{2(2x+4)-x-5}{4(-y+2x+4)}\right)}{2} =$\\ $log(x + 1)+\%c]$$

(%i13)	method;

$$abel2\leqno{(\%o13) }$$

(%i14)	eqn1:'diff(y,x)=(1-3*x-3*y)/(1+x+y);

$$\frac{d}{d\,x}\,y=\frac{-3\,y-3\,x+1}{y+x+1}\leqno{(\%o14) }$$

(%i15)	contrib_ode(eqn1,y,x);

$$[\frac{2\,log\left( y+x-1\right) +y+3\,x}{2}=\%c]\leqno{(\%o15) }$$

(%i16)	method;

$$lie\leqno{(\%o16) }$$

3.8.5.4 Решение однородных линейных уравнений

Другие полезные функции пакета contrib_ode: $$odelin$$ и $$ode_check$$.

Функция $$odelin$$ решает однородные линейные уравнения первого и второго порядка, и возвращает фундаментальное решение ОДУ.

Пример:

(%i4)	odelin(x*(x+1)*'diff(y,x,2)+(x+5)*'diff(y,x,1)+(-4)*y,y,x);

...trying factor method...
solving 7 equations in 4 variables...
trying the Bessel solver...solving 1 equations in 2 variables...
trying the F01 solver...
solving 1 equations in 3 variables...
trying the spherodial wave solver...
solving 1 equations in 4 variables...
trying the square root Bessel solver...
solving 1 equations in 2 variables...
trying the 2F1 solver...
solving 9 equations in 5 variables

$${\frac{gauss\_a\left( -6,-2,-3,-x\right) }{{x}^{4}},\frac{gauss\_b\left( -6,-2,-3,-x\right) }{{x}^{4}}}\leqno{(\%o4) }$$

Примечание: функции $$gauss_a$$ и $$gauss_b$$ — специальные функции, представляющие собой решения гипергеометрического уравнения.

Функция $$ode_check$$ позволяет подставить в ОДУ найденное решение.

Пример:

(%i1)	load("contrib_ode")$
(%i2)	eqn:(1+x^2)*'diff(y,x,2)-2*x*'diff(y,x);

$$\left( {x}^{2}+1\right) \,\left( \frac{{d}^{2}}{d\,{x}^{2}}\,y\right) -2\,x\,\left( \frac{d}{d\,x}\,y\right) \leqno{(\%o2) }$$

(%i3)	odelin(eqn,y,x);
...trying factor method...solving 7 equations in 4 variables

$${1,x\,\left( {x}^{2}+3\right) }\leqno{(\%o3) }$$

(%i4)	ode_check(eqn,y=x*(x^2+3));

$$0\leqno{(\%o4) }$$

(%i5)	ode_check(eqn,y=1);

$$0\leqno{(\%o5) }$$

3.8.6 Численные методы решения ОДУ

Однако в ряде случаев отыскать символьное решение ОДУ в достаточно компактном виде невозможно. В этом случае целесообразно использовать численные методы. Maxima включает пакет расширения dynamics, позволяющий проинтегрировать систему ОДУ методом Рунге-Кутта.

Начиная с версии 5.12, Maxima включает пакет dynamics (его необходимо загружать перед использованием). Помимо метода Рунге- Кутта, пакет dynamics включает ряд функций для построения различных фракталов.

Метод Рунге-Кутта реализует функция $$rk$$. Синтаксис вызова её вызова: $$rk([eq], [vars], [init], [t_{range}])$$, где $$eq$$ — список правых частей уравнений; $$vars$$ — список зависимых переменных; $$init$$ — список начальных значений; $$t_{range}$$ — список $$[t,t_0,t_{end},ht]$$, содержащий символьное обозначение независимой переменной $$(t)$$, её начальное значение $$(t_0)$$, конечное значение $$(t_{end})$$, шаг интегрирования $$(ht)$$.

Пример:

Решить ОДУ $$\frac{dx}{dt} = 4x^2-4y^2; \qquad \frac{dy}{dt} = y^2-x^2+1;$$ при $$t = [0\dots4],\ x(0)=-1,25,\ y(0)=0,75$$.

Используем пакет dynamics.

(%i1)	load(''dynamics'')$

Выбираем шаг интегрирования 0,02.

(%i2)	sol:rk([4*x^2-4*y^2,y^2-x^2+1],[x,y],
	[-1.25,0.75],[t,0,4,0.02]);

В результате решения получаем список значений в формате [[$$t,x,y$$]].

(%i1)	load("dynamics")$
(%i2)	rp1:4*x^2-4*y^2;

$$4\,{x}^{2}-4\,{y}^{2}\leqno{(\%o2) }$$

(%i3)	rp2:y^2-x^2+1;

(рис 3.16) Пример графического решения системы ОДУ численным методом $${y}^{2}-{x}^{2}+1\leqno{(\%o3) }$$

(%i4)	sol:rk([rp1,rp2],[x,y],[-1.25,0.75],[t,0,4,0.02])$

Список $$sol$$ не выводим на экран (он достаточно длинный, поэтому завершаем ввод команды символом $).

Для построения графика решения преобразуем полученный список, построив отдельно список значений $$t$$(список $$xg$$ в примере), $$x$$(список $$yg1$$), $$y$$(список $$yg2$$). При построении графика используем опцию $$discrete$$.

(%i5)	len:length(sol);

$$201\leqno{(\%o5) }$$

(%i6)	xg:makelist(sol[k][1],k,1,len)$
(%i7)	yg1:makelist(sol[k][2],k,1,len)$
(%i8)	yg2:makelist(sol[k][3],k,1,len)$
(%i9)	plot2d([[discrete,xg,yg1],[discrete,xg,yg2]]);

Результат решения представлен на рис. 3.16

Аналогичный, хотя и несколько более сложный пример — моделирование аттрактора Лоренца.

3.9 Ряды Фурье по ортогональным системам

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

3.9.1 Понятие ряда Фурье

Пусть даны две функции $$f(x)$$ и $$g(x)$$, произведение которых интегрируемо на отрезке [$$a,b$$]. Функции $$f(x)$$ и $$g(x)$$, называются ортогональными на [$$a,b$$], если выполняется условие $$\int _{a}^{b}f(x)g(x) \rho(x)dx=0,$$ где $$\rho(x)$$ — весовая функция.

Функциональная последовательность $$\{\varphi _{n} (x)\}=\{\varphi _{0} (x),\varphi _{1} (x),\dots\varphi _{n} (x),\dots\}$$ называется ортогональной на [$$a,b$$], если выполняется условие: $$\int _{a}^{b}\varphi _{n} \left(x\right) \varphi _{m} \left(x\right) \rho \left(x\right) dx=0, \forall n\ne m.$$.

Функциональная последовательность $$\left\{\varphi _{n} \left(x\right)\right\}$$ называется ортонормированной на [$$a,b$$], если $$\int_{a}^{b}\varphi_{n} \left(x\right)\varphi_{m} \left(x \right) \rho\left(x\right)dx = \begin{cases} 1, \mbox{если } \, n=m \\ 0, \mbox{если } \, n \neq m \end{cases}$$

Часто используемая последовательность из тригонометрических функций $$1,\cos(x),\sin(x),\cos(2x),\sin(2x),... ,\cos(nx),\sin(nx),..$$. ортогональна на отрезке $$[-\pi ,\pi]$$ с весовой функцией $$\rho(x)=1$$.

Проверим свойство ортогональности, вычисляя соответствующие интегралы. При $$m \ne n$$ получаем:

$$\begin{align*} \int _{-\pi }^{\pi }1\cdot \sin (nx)dx = -\left. \frac{\cos (nx)}{n} \right|_{-\pi }^{\pi } =0 ,\\ \int _{-\pi }^{\pi }1\cdot \sin (nx)dx = - \left. \frac{\cos (nx)}{n} \right|_{-\pi }^{\pi } =0 ,\\ \int _{-\pi }^{\pi }\cos (nx)dx = \left. \frac{1}{n} \sin (nx)\right|_{-\pi }^{\pi } =0, \forall n\in N; \\ \int _{-\pi }^{\pi }\cos (mx)\cdot \cos (nx)dx = \\ \frac{1}{2} \int _{-\pi }^{\pi }(\cos ((m-n)x)+ \cos ((m+n)x))dx= \\ \left. \frac{1}{2} \left(\frac{\sin ((m-n)x)}{m-n} +\frac{\sin (m+n)x}{m+n}\right) \right|_{-\pi }^{\pi } =0 \end{align*}$$

Если же $$m = n$$, то

$$\begin{align*} \int _{-\pi }^{\pi }\cos ^{2} (mx)dx=\frac{1}{2} \int _{-\pi }^{\pi }\left(1+\cos (2mx)\right)dx= \\ \frac{1}{2} \left. \left(x+\frac{\sin (2mx)}{2m} \right)\right|_{-\pi }^{\pi } =\pi \end{align*}$$

Следовательно, $$\int _{-\pi }^{\pi }\cos (mx)\cos (nx)dx=\left\{_{\pi ,\ m=n.}^{0,\ m\ne n;}\right$$.. Аналогичным образом устанавливаем, что $$\int _{-\pi }^{\pi }\sin (mx)\sin (nx)dx=\left\{_{\pi ,\ m=n.}^{0,\ m\ne n;}\right$$..

Остаётся вычислить интеграл $$\int _{-\pi }^{\pi }\cos (mx)\sin (nx)dx$$..

Поскольку подинтегральная функция является нечётной, то

$$\int _{-\pi }^{\pi }\cos (mx)\sin (nx)dx=0,$$

Как следует из приведённых равенств, любые две различные функции тригонометрической последовательности ортогональны на отрезке $$[-\pi ,\pi]$$.

Другой широко используемой последовательностью ортогональный функций является последовательность полиномов Лежандра. Полином Лежандра степени n можно представить через формулу Родрига в виде:

$$P_n(z)=\frac{1}{2^n n!}\frac{d^n}{dz^n}(z^2-1)^n .$$

Они также могут быть вычислены по рекуррентной формуле:

$$P_{n+1}(x)=\frac{2n+1}{n+1}xP_n(x)-\frac{n}{n+1}P_{n-1}(x).$$

Полиномы Лежандра ортогональны на отрезке [-1, 1] с весом $$\rho(x)=1$$:

$$\int_{-1}^{1}P_k(x)P_l(x)\,dx = \begin{cases} \frac{2}{2k+1}, \mbox{если } \, k=l \\ 0, \mbox{если } \, k \neq l \end{cases}.$$

Ещё одной важной последовательностью ортогональных функций является последовательность полиномов Чебышёва. Полиномы Чебышёва первого рода $$T_n(x)$$ степени n можно определить с помощью равенства:

$$T_n(\cos(\theta))=\cos(n\theta),$$

или, что почти эквивалентно,

$$T_n(z)=\cos(n \arccos(z)).$$

Они также могут быть вычислены по рекуррентной формуле:

$$\begin{align*} T_0(x)=1\\ T_1(x)=x\\ T_{n+1}(x) = 2xT_n(x) - T_{n-1}(x). \end{align*}$$

Полиномы Чебышёва ортогональны на отрезке [-1, 1] с весом $$\rho(x)=\cfrac{1}{\sqrt{1-x^2}}$$:

$$\int_{-1}^{1}T_k(x)T_l(x) \cfrac{1}{\sqrt{1-x^2}} \,dx = \begin{cases} \cfrac{\pi}{2}, \mbox{если } \, k=l\neq 0 \\ \pi, \mbox{если } \, k=l=0 \\ 0, \mbox{если } \, k \neq l \end{cases}$$

3.9.2 Вычисление коэффициентов тригонометрических рядов Фурье

Члены тригонометрического ряда $$\cfrac{a_{0} }{2} +\displaystyle\sum _{n=1}^{\infty }\left(a_{n} \cos (nx)+b_{n} \sin (nx)\right)$$ являются периодическими функциями с общим периодом $$2\pi$$, поэтому и сумма этого ряда $$S(x)$$ также будет периодической функцией с периодом $$2\pi$$.

Предположим, что $$2\pi$$–периодическую функцию $$f(x)$$ можно разложить в тригонометрический ряд, равномерно сходящийся на отрезке $$[-\pi, \pi]$$.

$$\left(x\right)=\frac{a_{0} }{2} +\sum _{n=1}^{\infty }\left(a_{n} \cos (nx)+b_{n} \sin (nx)\right)$$

Рассмотрим вопрос об определении коэффициентов $$a_0, a_n$$ и $$b_n (n =1, 2,...)$$. Для этого применим теорему о почленном интегрировании функционального ряда. Проинтегрируем обе части равенства в пределах от $$-\pi$$ до $$\pi$$:

$$\[\int _{-\pi }^{\pi }f\left(x\right)dx=\frac{a_{0} }{2} \int _{-\pi }^{\pi }dx+\sum _{n=1}^{\infty }\left(a_{n} \int _{-\pi }^{\pi }\cos (nx)dx+b_{n} \int _{-\pi }^{\pi }\sin (nx)dx \right) .\]$$

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

$$\[\int _{-\pi }^{\pi }f(x)dx=\pi a_{0} .\]$$

Следовательно,

$$a_{0} =\frac{1}{\pi } \int _{-\pi }^{\pi }f\left(x\right)dx .$$

Для того чтобы найти $$a_n (n = 1, 2,...)$$, обе части этого равенства умножим на cos(mx) и проинтегрируем на отрезке $$[-\pi,\pi]$$. Поскольку система тригонометрических функций ортогональна, то $$\[\int _{-\pi }^{\pi }\cos (mx)\cos (nx)dx=0,\ \ \int _{-\pi }^{\pi }\cos (mx)\sin (nx)dx=0 \]$$ для $$\forall m,n\in {\mathbb N}$$, если $$m \ne n$$.

Это означает что все с интегралы, встречающиеся в правой части, будут равны нулю, исключение составляет интеграл, который получается при $$m = n$$. Этот интеграл равен $$\pi$$. Поэтому $$\[\int _{-\pi }^{\pi }f\left(x\right)\cos (nx)dx=a_{n} \int _{-\pi }^{\pi }\cos ^{2} (nx)dx=\pi a_{n} ,\]$$ откуда $$a_{n} =\cfrac{1}{\pi } \displaystyle\int _{-\pi }^{\pi }f\left(x\right)\cos (nx)dx, \, n =1,2,\dots$$

Аналогично, умножив обе части равенства на $$\sin(mx)$$ и проинтегрировав на отрезке $$[-\pi; \pi]$$, получаем, что $$b_{n} =\cfrac{1}{\pi } \displaystyle\int _{-\pi }^{\pi }f\left(x\right)\sin ( nx)dx,\ n=1,2,\dots$$

Итак, если функцию $$f(x)$$ можно представить в виде тригонометрического ряда, то коэффициенты $$a_0, a_n, b_n$$ вычисляются по приведённым формулам и называются коэффициентами Фурье для функции $$f(x)$$ (а ряд — соответственно рядом Фурье для $$f(x)$$).

Промежуток интегрирования $$[-\pi,\pi]$$ для периодической с периодом $$2\pi$$ функции можно заменить любым промежутком $$[a,a + 2\pi]$$, $$a\in \mathbb R$$, длина которого равна $$2\pi$$.

Функция $$f(x)$$ называется кусочно-гладкой на отрезке [$$a,b$$] если функция $$f(x)$$ и её производная на [$$a,b$$] имеют конечное число точек разрыва первого рода.

Достаточные условия разложимости функции в ряд Фурье даёт теорема Дирихле: если $$f(x)$$ — периодическая с периодом $$2\pi$$ кусочногладкая на $$[-\pi; \pi]$$ функция, то её ряд Фурье сходится в любой точке этого отрезка и его сумма равна:

значению функции $$f(x)$$, когда $$x$$ — точка непрерывности функции $$f(x)$$;

$$\cfrac{f\left(x-0\right)+f(x+0)}{2}$$, когда $$x$$ — точка разрыва функции $$f(x)$$, при этом $$\cfrac{f(x-0)+f(x+0)}{2} =\cfrac{a_{0} }{2} +\displaystyle\sum _{n=1}^{\infty }\left(a_{n}\cos (nx)+b_{n} \sin (nx)\right) .$$

Отметим, что на практике чаще всего встречаются функции, которые удовлетворяют условиям теоремы Дирихле.

Пример: периодическую с периодом $$2\pi$$ функцию $$f(x) = x$, $ -\pi < x < \pi$$ разложить в ряд Фурье.

Вычислим коэффициенты Фурье (используем Maxima):

(%i1)	n:5;

$$5\leqno{(\%o1) }$$

(%i2)	f(x):=x;

$$f\left( x\right) :=x\leqno{(\%o2) }$$

(%i3)	a0:1/%pi*integrate(f(x),x,-%pi,%pi);

$$0\leqno{(\%o3) }$$

(%i4)
for k:1 thru n do a[k]:1/%pi*integrate(f(x)*cos(k*x),x,-%pi,%pi);

$$done\leqno{(\%o4) }$$

(%i5)
for k:1 thru n do b[k]:1/%pi*integrate(f(x)*sin(k*x),x,-%pi,%pi);

$$done\leqno{(\%o5) }$$

(%i6)	for k:1 thru n do display(a[k],b[k]);

$${a}_{1}=0\ {b}_{1}=2\ {a}_{2}=0\ {b}_{2}=-1\ {a}_{3}=0\ {b}_{3}=\frac{2}{3}\ {a}_{4}=0\ {b}_{4}=-\frac{1}{2}\ {a}_{5}=0\ {b}_{5}=\frac{2}{5}$$ $$done\leqno{(\%o6) }$$

(%i7)
fun(x):=a0/2+sum(a[k]*cos(k*x),k,1,n)+sum(b[k]*sin(k*x),k,1,n);

$$fun\left( x\right) :=\frac{a_0}{2}+sum\left( {a}_{k}\,cos\left( k\,x\right) ,k,1,n\right) +sum\left( {b}_{k}\,sin\left( k\,x\right) ,k,1,n\right) \leqno{(\%o7) }$$

(%i8)	wxplot2d([f(x),fun(x)], [x,-5,5], [nticks,20]);

Данная функция $$f(x)$$ удовлетворяет условиям теоремы Дирихле, её график в сравнении с графиком частичной суммы ряда Фурье $$fun(x)$$ изображён на рис. 3.17.

3.9.3 Ряды Фурье для чётных и нечётных функций

Предположим, что $$f(x)$$ — нечётная $$2\pi$$–периодическая функция. В этом случае $$f(x)\cos(nx)$$ — чётная функция, поскольку верно равенство$$ f(-x)\cos(-nx) = f(x)\cos(nx)$$, a $$f(x)\sin(nx)$$ — нечётная функция, так как $$(-x)\sin(-nx) = -f(x)\sin(nx)$$ Поэтому коэффициент ряда Фурье $$a_n, b_n$$ равны:

$$a_{n} =\cfrac{1}{\pi } \displaystyle\int_{-\pi }^{\pi }f(x)\cos (nx)dx=\cfrac{2}{\pi } \int _{0}^{\pi }f(x)\cos (nx)dx\ (n=0,1,\dots ),$$ $$b_{n} =\cfrac{1}{\pi } \displaystyle\int _{-\pi }^{\pi }f(x)\sin (nx)dx=0\ (n=1,2,\dots ).$$ (рис 3.17) График функции y = f(x) и суммы первых пяти членов ряда Фурье

Следовательно, ряд Фурье чётной функции содержит только косинусы, т.е. $$f(x)=\cfrac{a_{0} }{2} +\displaystyle\sum _{n=1}^{\infty }a_{n} \cos (nx)$$.. Аналогично, если $$f(x)$$ — нечётная функция, то $$f(x)\cos(nx)$$ — нечётная, а $$f(x)\sin(nx)$$ — чётная функция.

Поэтому $$a_{n} =\cfrac{1}{\pi } \displaystyle\int _{-\pi }^{\pi }f(x)\cos (nx)=0\ (n=0,1,\dots ),\\ b_{n} =\cfrac{1}{\pi } \displaystyle\int_{-\pi }^{\pi }f(x)\sin (nx)dx= \cfrac{2}{\pi } \displaystyle\int_{0}^{\pi }f(x)\sin (nx)dx\ (n=1,2,\dots )$$

Следовательно, ряд Фурье нечётной функции содержит только синусы, т.е. $$f(x)=\displaystyle\sum _{n=1}^{\infty }b_{n} \sin (nx)$$.

Пример: Разложить в ряд Фурье периодическую с периодом $$2\pi$$ функцию, заданную на отрезке $$[-\pi,\pi]$$ равенством $$f(x) = x^2$$.

(рис 3.18) График функции (точки) и суммы первых пяти членов ряда Фурье (сплошная линия)

Данная функция $$y=x^2$$ является чётной (рис. 3.18), поэтому её ряд Фурье содержит только косинусы. Вычисляем коэффициенты этого ряда: $$b_n = 0, n = 1, 2,...$$

Для вычисления коэффициентов $$a_n$$ ряда Фурье создаём функцию $$fun$$, входными параметрами которой являются имя независимой переменной (в примере это $$x$$), число суммируемых членов ряда ($$n$$, в дальнейшем функция вызывается при $$n = 5$$) и символьное выражение, определяющее функцию, для которой строится разложение ($$f$$, функция $$fun$$ вызывается с $$f = x^2$$).

$$Пример:$$

(%i1)	fun(x,n,f):=(for k:0 thru n do
	a[k]:1/%pi*integrate(f*cos(k*x),x,-%pi,%pi),
	a[0]/2 +sum(a[k]*cos(k*x),k,1,n))$

(%i2)	fun(x,5,x^2);

$$(\%o2) -\frac{4\,\mathrm{cos}\left( 5\,x\right) }{25}+\frac{\mathrm{cos}\left( 4\,x\right) }{4}-\frac{4\,\mathrm{cos}\left( 3\,x\right) }{9}+\mathrm{cos}\left( 2\,x\right) -4\,\mathrm{cos}\left( x\right) +\frac{{\pi }^{2}}{3}$$

Для аналитического вычисления коэффициентов ряда Фурье функции $$y = |x|$$ функцию $$fun$$ необходимо немного изменить, предусмотрев различные выражения для подинтегрального выражения на полуинтервалах $$[-\pi, 0) и (0,\pi]$$ (выражения $$f_1$$ и $$f_2$$ в списке параметров функции). Текст программы на макроязыке Maxima:

fun12(x,n,f1,f2):=(for k:0 thru n do
a[k]:1/%pi*(integrate(f1*cos(k*x),x,-%pi,0)+
		integrate(f2*cos(k*x),x,0,%pi)),
a[0]/2+sum(a[k]*cos(k*x),k,1,n))$

(рис 3.19) График функции y = |x| (точки) и суммы первых пяти членов ряда Фурье (сплошная линия)

Функция является $$y = |x|$$ также является чётной (рис. 3.19), поэтому её ряд Фурье содержит только косинусы.

Результаты вычисления коэффициентов ряда Фурье для этой функции:

(%i1)	fun12(x,5,-x,x);

$$-\frac{4\,\mathrm{cos}\left( 5\,x\right) }{25\,\pi }-\frac{4\,\mathrm{cos}\left( 3\,x\right) }{9\,\pi }-\frac{4\,\mathrm{cos}\left( x\right) }{\pi }+\frac{\pi }{2}\leqno{(\%o1)}$$

Для построения графика функции $$y = |x|$$ создаём функцию $$fg(x)$$, которая использована для построения графика на рис. 3.19.

(%i3)	fg(x):=if x>0 then x else -x$

3.9.4 Разложение функций в ряд Фурье на отрезке [0,π]

Пусть $$f(x)$$ определена на отрезке $$[0,\pi]$$. Для того, чтобы функцию $$f(x)$$ разложить в ряд Фурье на этом отрезке, доопределим эту функцию произвольным образом на интервале $$[-\pi, 0]$$. Рассмотрим два случая:

Функцию $$f(x)$$, заданную на $$[0,\pi]$$, продолжим на интервал $$[-\pi, 0]$$ так, что вновь полученная функция $$f_1(x)$$ , была чётной:

$$f_{1} = \begin{cases} f(-x), \mbox{если } \, x\in [-\pi ,0] \\ f(x), \mbox{если } \, x\in [0,\pi ] \end{cases}.$$

В таком случае говорят, что $$f(x)$$ продолжена на $$[-\pi, 0]$$ чётным образом. Поскольку $$f_1(x)$$ — чётная на $$[-\pi,\pi]$$ функция, то её ряд Фурье содержит только косинусы:

$$\[f_{1} (x)=\frac{a_{0} }{2} +\sum _{n=1}^{\infty }a_{n} \cos (nx).\]$$

Поскольку на отрезке $$[0,\pi]$$ имеет место равенство $$f_1(x) = f(x)$$, то ряд Фурье для функции $$f_1(x)$$ будет и рядом Фурье для $$f(x)$$ на $$[0,\pi]$$

Функцию $$f(x)$$, заданную на $$[0,\pi]$$, продолжим на интервал $$[-\pi, 0]$$ нечётным образом:

$$f_{2} = \begin{cases} -f(-x), \mbox{если } \, x\in [-\pi ,0[ \\ f(x), \mbox{если } \, x\in [0,\pi] \end{cases}.$$

Поскольку $$f_2(x)$$ — нечётная на $$[-\pi,\pi]$$ функция, то её ряд Фурье содержит только синусы:

$$\[f_{2} (x)=\sum _{n=1}^{\infty }b_{n} \sin (nx).\]$$

Так как $$f_2(x) = f(x)$$ при $$\forall x\in [0,\pi]$$, то полученный ряд Фурье для $$f_2(x)$$ и будет рядом Фурье для $$f(x)$$ на $$[0,\pi]$$.

Пример: Функцию $$f(x) = 2x + 1$$, определённую на отрезке $$[0,\pi]$$, разложить в ряд Фурье: 1)по косинусам; 2)по синусам.

1) Функцию $$f(x)$$ продолжим на $$[-\pi, 0]$$ чётным образом, т.е. составим новую функцию $$f_1(x)$$ по формуле:

$$f_{1}(x) = \begin{cases} -2x+1, \mbox{если } \, x\in [-\pi ,0[ \\ 2x+1, \mbox{если } \, x\in [0,\pi] \end{cases}.$$

Вычисляем коэффициенты Фурье для этой функции при помощи функции $$fun12$$:

(%i1)	fleft:-2*x+1;

$$(\%o1)\quad 1-2\,x$$

(%i2	 fright:2*x+1;

$$(\%o2)\quad 2\,x+1$$

(%i3)	funcos(x,7,fleft,fright);

$$(\%o3)\ -\frac{8\,cos\left( 7\,x\right) }{49\,\pi }-\frac{8\,cos\left( 5\,x\right) }{25\,\pi }-\frac{8\,cos\left( 3\,x\right) }{9\,\pi }-\frac{8\,cos\left( x\right) }{\pi }+\frac{2\,{\pi }^{2}+2\,\pi }{2\,\pi }$$ (рис 3.20) График функции y = 2x+1, продолженной чётным образом, и суммы семи членов соответствующего ряда

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

2) Функцию $$f(x)$$ продолжим на $$[-\pi, 0]$$ нечётным образом. Составим новую функцию $$f_2(x)$$ по формуле $$f_{2} (x)=\left\{_{2x+1,\ x\in [0,\pi]}^{2x-1,\ x\in [-\pi ,0[} \right.$$

Вычислим коэффициенты Фурье для этой функции, используя функцию $$fun12sin$$, аналогичную приведённой выше.

Пример:

(%i1)	fleft:2*x-1$
(%i2)	fright:2*x+1$
(%i3)	f(x):=(if x>0 then fright else fleft)$
(%i4)	fun12sin(x,n,f1,f2):=(for k:1 thru n do
	b[k]:1/%pi*(integrate(f1*sin(k*x),x,-%pi,0)
	+integrate(f2*sin(k*x),x,0,%pi)),
	sum(b[k]*sin(k*x),k,1,n))$
(%i5)	fun12sin(x,7,fleft,fright);

$$(\%o5)\ \frac{\left( \frac{2\,\left( 2\,\pi +1\right) }{7}+\frac{2}{7}\right) \,sin\left( 7\,x\right) }{\pi }+\frac{\left( \frac{1}{3}-\frac{2\,\pi +1}{3}\right) \,sin\left( 6\,x\right) }{\pi }+ \frac{\left( \frac{2\,\left( 2\,\pi +1\right) }{5}+\frac{2}{5}\right) \,sin\left( 5\,x\right) }{\pi }+\\ \frac{\left( \frac{1}{2}-\frac{2\,\pi +1}{2}\right) \,sin\left( 4\,x\right) }{\pi }+\frac{\left( \frac{2\,\left( 2\,\pi +1\right) }{3}+\frac{2}{3}\right) \,sin\left( 3\,x\right) }{\pi }-2\,sin\left( 2\,x\right) +\frac{\left( 4\,\pi +4\right) \,sin\left( x\right) }{\pi }$$

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

(рис 3.21) Сравнение графика функции y = 2x + 1 при нечётном продолжении и суммы семи членов соответствующего ряда Фурье

3.9.5 Ряд Фурье для функций с периодом 2l

Пусть $$f(x)$$ — периодичная с периодом $$2\ell\ \left(\ell \ne \pi \right)$$ функция, которая на отрезке $$[-\ell ,\ell]$$ удовлетворяет условиям теоремы Дирихле. Разложим её на этом отрезке в ряд Фурье. Обозначим

$$x=\frac{\ell t}{\pi } .$$

Тогда

$$\[f(x)=f\left(\frac{\ell t}{\pi } \right)=\varphi (t)\]$$

Функция $$\varphi (t)$$ — уже $$2\pi$$–периодическая функция, так как

$$\[\varphi (t+2\pi )=f\left(\frac{\ell }{\pi } (t+2\pi )\right)=f\left(\frac{\ell t}{\pi } +2\ell \right)=f\left(\frac{\ell t}{\pi } \right)=\varphi (t).\]$$

Функцию $$\varphi (t)$$ разложим в ряд Фурье на отрезке $$[-\pi ,\pi]$$

$$\varphi (t)=f\left(\frac{\ell t}{\pi } \right)=\frac{a_{0} }{2} +\sum _{n=1}^{\infty }\left(a_{n} \cos (nt)+b_{n} \sin (nt)\right) \dots$$

Коэффициенты этого ряда вычисляются по формулам:

$$a_{n} =\frac{1}{\pi } \int _{-\pi }^{\pi }f\left(\frac{\ell t}{\pi } \right) \cos (nt)dt,n=0,1,\dots ,$$ $$b_{n} =\frac{1}{\pi } \int _{-\pi }^{\pi }f\left(\frac{\ell t}{\pi } \right) \sin (nt)dt,n=1,2,\dots$$

Возвращаясь к прежней переменной $$x$$, из равенства (3.3) имеем $$t=\frac{\pi x}{\ell }$$. Тогда ряд (3.4) можно представить в виде

$$f(x)=\frac{a_{0} }{2} +\sum _{n=1}^{\infty }\left(a_{n} \cos (\frac{n\pi x}{\ell } )+b_{n} \sin (\frac{n\pi x}{\ell } )\right).$$

В интегралах (3.5) и (3.6) произведём замену переменной:

$$a_{n} =\frac{1}{\pi } \int _{-\pi }^{\pi }f\left(\frac{\ell t}{\pi } \right)\cos (nt)dt= \frac{1}{\ell } \int _{-\ell }^{\ell }f(x)\cos (\frac{n\pi x}{\ell } )dx,n=0,1,\dots\\ b_{n} =\cfrac{1}{\pi } \displaystyle\int _{-\pi }^{\pi }f\left(\frac{\ell t}{\pi } \right)\sin (t)dt=\cfrac{1}{\ell } \displaystyle\int _{-\ell }^{\ell }f(x)\sin (\frac{n\pi x}{\ell } )dx, n=1,2,\dots.$$

Если $$f(x)$$ — чётная на $$\left[-\ell ,\ell \right]$$ функция, то $$b_n = 0 (n = 1, 2,...)$$, а $$a_{n} = \cfrac{2}{\ell } \displaystyle\int _{0}^{\ell }f(x)\cos (\frac{n\pi x}{\ell } )dx,\ (n=0,1,\dots),$$, ряд Фурье такой функции имеет вид:

$$\[f(x)=\frac{a_{0} }{2} +\sum _{n=1}^{\infty }a_{n} \cos (\frac{n\pi x}{\ell } ) .\]$$

Если $$f(x)$$ — нечётная на $$[-\ell ,\ell ]$$ функция, то $$a_n = 0 (n = 0, 1, 2,...)$$, $$b_{n} =\cfrac{2}{\ell } \displaystyle\int _{0}^{\ell }f(x)\sin (\frac{n\pi x}{\ell } ) dx,\ (n=1,2,\dots )$$, ряд Фурье имеет вид:

$$\[f(x)=\sum _{n=1}^{\infty } b_{n} \sin (\frac{n\pi x}{\ell } )dx.\]$$

Пример: Разложить в ряд Фурье периодическую с периодом $$T =2$$ функцию $$f(x)$$, заданную формулой

$$f(x) = \begin{cases} x, \mbox{если } \, 0<x\le 1 \\ 0, \mbox{если } \, -1<x\le 0 \end{cases}.$$

Эта функция на отрезке [-1, 1] удовлетворяет условиям теоремы Дирихле. Ряд Фурье для данной функции:

$$\[f(x)=\frac{1}{4} -\frac{2}{\pi ^{2} } \sum _{k=0}^{\infty }\frac{\cos ((2k+1)x)}{(2k+1)^{2} } +\frac{1}{\pi } \sum _{k=1}^{\infty }\frac{(-1)^{k+1} }{k} \sin (k\pi x)\]$$

Сумма этого ряда в точках $$x=\pm 1,\pm 3,\dots$$ равна $$\cfrac{1}{2}$$.

Рассмотрим видоизменение функции Maxima, необходимой для вычисления коэффициентов ряда Фурье для функции с периодом $$[-\ell, \ell]$$. Рассмотрим текст функции $$fun12l$$:

fun12l(x,n,l,f1,f2):=(for k:0 thru n do
	a[k]:1/l*(integrate(f1*cos(%pi*k*x/l),x,-l,0)
	+integrate(f2*cos(%pi*k*x/l),x,0,l)),
for k:1 thru n do b[k]:1/l*(integrate(f1*sin(%pi*k*x/l),x,-l,0)+
	integrate(f2*sin(%pi*k*x/l),x,0,l)),
	a[0]/2+sum(a[k]*cos(%pi*k*x/l),k,1,n)+
	sum(b[k]*sin(%pi*k*x/l),k,1,n))$

Основное изменение по сравнение с вариантами, приведёнными выше — использование тригонометрических функций $$\sin \left( \frac{\pi k x}{\ell} \right)$$ и $$\sin \left( \frac{\pi k x}{\ell} \right)$$

Вывод Maxima для первых семи членов ряда Фурье:

(%i6)	fun12l(x,7,1,0,x);

$$(\%o6)\ \frac{sin\left( 7\,\pi \,x\right) }{7\,\pi }-\frac{2\,cos\left( 7\,\pi \,x\right) }{49\,{\pi }^{2}}-\frac{sin\left( 6\,\pi \,x\right) }{6\,\pi }+\frac{sin\left( 5\,\pi \,x\right) }{5\,\pi }-\frac{2\,cos\left( 5\,\pi \,x\right) }{25\,{\pi }^{2}}-\frac{sin\left( 4\,\pi \,x\right) }{4\,\pi }+\frac{sin\left( 3\,\pi \,x\right) }{3\,\pi }-\frac{2\,cos\left( 3\,\pi \,x\right) }{9\,{\pi }^{2}}-\frac{sin\left( 2\,\pi \,x\right) }{2\,\pi }+\frac{sin\left( \pi \,x\right) }{\pi }-\frac{2\,cos\left( \pi \,x\right) }{{\pi }^{2}}+\frac{1}{4}$$

Для построения графика собственно анализируемой функции (её представляет кусочно-непрерывная функция $$f(x)$$) и частичной суммы её ряда Фурье из результатов разложения формируем новую функцию $$g(x)$$, после чего стандартной командой строим график:

(%i7)	g(x):="%$
(%i8)	f(x):=(if x<0 then 0 else x)$
(%i9)	wxplot2d([g(x),f(x)], [x,-2.2,1.6]);

Графическая иллюстрация, показывающая сопоставление рассматриваемой функции и ряда Фурье на заданном отрезке — на рис. 3.22.

(рис 3.22) График функции f(x) и суммы первых семи членов ряда Фурье

3.9.6 Комплексная форма ряда Фурье

Пусть функция $$f(x)$$ на $$\left[-\pi ,\pi \right]$$ разложена в ряд Фурье

$$f(x)=\frac{a_{0} }{2} +\sum _{n=1}^{\infty }\left(a_{n} \cos (x)+b_{n} \sin (x)\right) .$$

Воспользуемся формулами Эйлера:

$$\[\cos (nx)=\frac{e^{inx} +e^{-inx} }{2} ,\ \sin (nx)=\frac{e^{inx} -e^{-inx} }{2i} .\]$$

Подставим эти выражения в ряд (3.8), имеем:

$$f(x)=\frac{a_{0} }{2} +\sum _{n=1}^{\infty} \left(a_{n} \frac{e^{inx} +e^{-inx} }{2} + b_{n} \frac{e^{inx} -e^{-inx} }{2i} \right) = \\ =\frac{a_{0} }{2} + \sum _{n=1}^{\infty }\left( a_{n} \frac{e^{inx}+e^{-inx} }{2} - ib_{n} \frac{e^{inx} -e^{-inx} }{2} \right) = \\ =\frac{a_{0} }{2}+\sum _{n=1}^{\infty }\left( \frac{a_{n} -ib_{n} }{2} \cdot e^{inx} +\frac{a_{n} +ib_{n} }{2} \cdot e^{-inx} \right).$$

Обозначим:

$$\frac{a_{0} }{2} =c_{0} ,\ \frac{a_{n} -ib_{n} }{2} =c_{n} ,\ \frac{a_{n} +ib_{n} }{2} =c_{-n} .$$

Тогда

$$f(x) =c_{0} +\sum _{n=1}^{\infty }\left(c_{n} \cdot e^{inx} +c_{-n} e^{-inx} \right) = \\ = c_{0} +\sum _{n=1}^{\infty }c_{n} e^{inx} +\sum _{n=1}^{\infty }c_{-n} e^{-inx} = \\ = c_{0} +\sum _{n=1}^{\infty }c_{n} e^{inx} +\sum _{n=-\infty }^{n=-1}c_{n} e^{inx} =\sum _{n=-\infty }^{\infty }c_{n} e^{inx} .$$

Следовательно

$$f(x)=\sum _{n=-\infty }^{\infty }c_{n} e^{inx}$$

Выражение (3.10) называется комплексной формой ряда Фурье функции $$f(x)$$ с комплексными коэффициентами Фурье $$c_n$$. Коэффициенты Фурье $$c_n$$ вычисляются по формулам $$(n=0, \pm 1, \pm 2,\dots)$$:

$$c_{n} =\frac{1}{2} \left(a_{n} -ib_{n} \right)= \frac{1}{2\pi} \int _{-\pi }^{\pi }f(x)\left[\cos (nx)-i\sin (nx)\right]dx=\\ =\frac{1}{2\pi } \int _{-\pi }^{\pi }f(x)\left[\cos (-nx)+i\sin (-nx)\right]dx=\frac{1}{2\pi } \int _{-\pi }^{\pi }f(x)e^{-inx} dx.$$

Если $$f(x)$$ — периодическая с периодом $$2\ell$$ функция, то её комплексный ряд Фурье имеет вид: $$\[f(x)=\sum _{n=-\infty }^{\infty }c_{n} e^{\cfrac{in\pi x}{\ell } } ,\]$$ а коэффициенты Фурье определяются по формуле

$$\[c_{n} =\frac{1}{2\ell } \int _{-\ell }^{\ell }f(x )e^{\cfrac{-in\pi x}{\ell }}dx.\]$$

Пример: Разложить в ряд Фурье с комплексными коэффициентами периодическую с периодом $$l = 2$$ функцию, заданную на отрезке [-1, 1] равенством $$f(x) = x^2$$.

(%i1)	n:5$ f:x^2$ l:1$ c(k):=
	1/2/l*integrate(f*exp(-%i*%pi*k*x/l),x,-l,l)$
	z:makelist(k-6, k, 1, 2*n+1)$
	cr:makelist(c(z[k]),k,1,2*n+1)$
	fk:makelist(cr[k]*exp(%i*%pi*z[k]*x/l),k,1,2*n+1)$
	g:sum(fk[k],k,1,2*n+1)$
	gend:trigreduce(ratsimp(rectform(g)));

$$(\%o9) \\ \frac{-144\,\mathrm{cos}\left( 5\,\pi \,x\right) +225\,\mathrm{cos}\left( 4\,\pi \,x\right) -400\,\mathrm{cos}\left( 3\,\pi \,x\right) +900\,\mathrm{cos}\left( 2\,\pi \,x\right) -3600\,\mathrm{cos}\left( \pi \,x\right) +300\,{\pi }^{2}}{900\,{\pi }^{2}}$$

В данном примере члены частичной суммы ряда Фурье представляются списком. В представленном вычислении $$z = -5, -4,... , 4, 5$$. Список $$cr$$ содержит коэффициенты ряда в комплексной форме (при суммировании от $$-n$$ до $$n$$ индекс элемента ряда содержится в $$z[k]$$). Собственно члены ряда Фурье скомпонованы в список $$fk$$, после суммирования которого получаем сумму ряда (выражение $$g$$). Для построения графика $$g(x)$$ необходимо упростить выражение $$g$$ (см. пример, результат упрощения — выражение $$gend$$). Очевидно, что для просмотра промежуточных результатов (они довольно объёмные) терминальные символы $ можно заменить на ;.

3.9.7 Дополнительные возможности: пакет fourie

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

Функция $$fourier$$ позволяет вычислить коэффициенты ряда Фурье (синтаксис вызова: $$fourier(f,x,p)$$), которая возвращает список коэффициентов Фурье $$f(x)$$, определённых на интервале $$[-p,p]$$. Собственно ряд Фурье позволяет построить функция $$fourexpand$$ (синтаксис вызова $$fourexpand(l,x,p,limit)$$), которая конструирует и возвращает ряд Фурье, используя список коэффициентов Фурье $$l$$ ($$limit$$ может быть и бесконечным, равным $$inf$$).

Коэффициенты рядов Фурье по синусам и по косинусам вычисляются функциями $$fourcos(f,x,p) foursin(f,x,p)$$ (синтаксис и аналогичны функции $$fourier$$).

Вычисления и подстановка $$\cos n\pi$$ и $$\sin n\pi$$ осуществляется специальной функцией $$foursimp(l)$$. Управление подстановкой осуществляется посредством флагов $$sinnpiflag$$ и $$cosnpiflag$$ (если они установлены в $$true$$, вычисление и подстановка выполняются, это режим по умолчанию).

Для управления процессом разложения различных функций в ряд Фурье предусмотрены следующие функции:

$$remfun$$. Синтаксис вызова $$remfun(f,expr)$$ или $$remfun(f,expr,x)$$. Данная функция позволяет заменить все вхождения функции $$f(arg)$$ в выражении $$expr$$ на $$arg$$ (в форме $$remfun(f,expr,x)$$ замена осуществляется, только если $$arg$$ содержит $$x$$);

$$funp$$. Данная функция (синтаксис вызова $$funp(f,expr)$$ или $$funp(f,expr,x)$$) возвращает $$true$$, если выражение $$expr$$ содержит функцию $$f$$ или конкретно $$f(x)$$;

$$absint$$. Данная функция позволяет вычислить неопределённый или определённый интеграл абсолютных значений функции $$f$$ (её определение может включать выражения $$abs(x), abs(sin(x)), abs(a) * exp(-abs(b) * abs(x)$$). Синтаксис вызова $$absint(f,x,halfplane) (halfplane = (pos,neg,both)$$ — часть числовой оси), $$absint(f,x)$$ (неопределённый интеграл по положительной полуоси), $$absint(f,x,a,b)$$ (определённый интеграл).

Общую форму ряда Фурье (после подстановки и упрощения) позволяет построить функция $$totalfourier(f,x,p)$$.

Коэффициенты интеграла Фурье на интервале (- inf, inf) позволяет вычислить функция $$fourint(f,x)$$, интеграла по косинусам или синусам на интервале (0, inf) — функции $$fourintcos(f,x)$$ и $$fourintsin(f,x)$$ соответственно.

Для использования пакета fourie его необходимо предварительно загрузить командой load("fourie").

(рис 3.23) График частичной суммы ряда Фурье для функции f(x) = x, построенной при помощи пакета fourie

Примеры использования пакета fourie (график полученной функции приведён на рис. 3.23):

(%i1)	load("fourie")$ fourier(x,x,%pi);

$$(\%t2)\ {a}_{0}=0\\ (\%t3)\ {a}_{n}=0\\ (\%t4)\ {b}_{n}=\frac{2\,\left( \frac{\mathrm{sin}\left( \pi \,n\right) } {{n}^{2}}-\frac{\pi \,\mathrm{cos}\left( \pi \,n\right) }{n}\right) }{\pi }\\ (\%o4)\ [\%t2,\%t3,\%t4]$$

(%i5)	foursimp(%);

$$(\%t5)\ {a}_{0}=0\\ (\%t6)\ {a}_{n}=0\\ (\%t7)\ {b}_{n}=-\frac{2\,{\left( -1\right) }^{n}}{n}\\ (\%o7)\ [\%t5,\%t6,\%t7]$$

(%i8)	fourexpand(%,x,%pi,10);

$$(\%o8)\ -\frac{\mathrm{sin}\left( 10\,x\right) }{5}+\frac{2\,\mathrm{sin}\left( 9\,x\right) }{9}-\frac{\mathrm{sin}\left( 8\,x\right) }{4}+\frac{2\,\mathrm{sin}\left( 7\,x\right) }{7}-\frac{\mathrm{sin}\left( 6\,x\right) }{3}+\frac{2\,\mathrm{sin}\left( 5\,x\right) }{5}-\frac{\mathrm{sin}\left( 4\,x\right) }{2}+\frac{2\,\mathrm{sin}\left( 3\,x\right) }{3}-\mathrm{sin}\left( 2\,x\right) +2\,\mathrm{sin}\left( x\right)$$

3.9.8 Дополнительные возможности: обобщённые ряды Фурье

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

Вычисление значений ортогональных полиномов в Maxima осуществляется при помощи пакета orthopoly, который позволяет оперировать полиномами Чебышёва, Лежандра, Эрмита, Якоби и др., а также рядом сферических функций.

Интегрируемая на интервале (-1, 1) кусочно-непрерывная функция может быть представлена обобщённым рядом Фурье (в данном случае — по полиномам Лежандра): $$f(x)=\sum\limits_{n=0}^{\infty}c_{n}P_{n}(x),$$ где $$P_n(x)$$ — полином Лежандра степени $$n, c_n$$ — коэффициенты Фурье для разложения по полиномам Лежандра. Значения $$c_n$$ вычисляются по формуле: $$c_{n}=\frac{2n+1}{2}\int\limits_{-1}^{1}f(x)P_{n}(x)dx.$$

Пример вычисления разложения функции $$y = e^x$$ на интервале (-1, 1) в ряд по полиномам Лежандра представлен следующими командами:

(%i1)	load(orthopoly)$ n:5$ f:exp(x)$ l:1$
c(m):=(2*m+1)/2*integrate(f*legendre*p (m, x),x,-l,l)$
z:makelist(k-1, k, 1, n+1)$
cr:makelist(c(z[k]),k,1,n+1)$
fk:makelist(cr[k]*legendre_p (z[k], x),k,1,n+1)$
g:sum(fk[k],k,1,n+1)$

(рис 3.24) График частичной суммы обобщённого ряда Фурье для функции

График полученного выражения $$g$$ в сравнении с функцией $$e^x$$ показан на рис. 3.24.

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