Файл: Примеры решения математических задач в maple. Оглавление.pdf

1
Доля П.Г.
Харьковский Национальный Университет механико – математический факультет
2011 г.
Примеры решения математических задач в MAPLE.
Оглавление
1. Разложение функции в ряд Фурье ..................................................................... 1 2. Решение задачи Коши разложением в степенной ряд..................................... 4 3. Метод Галеркина решения краевых задач ОДУ ............................................ 12
Литература. ............................................................................................................ 16
1. Разложение функции в ряд Фурье
В Maple не существует специального пакета для построения разложения функций в ряды Фурье. Однако разработка процедуры, выполняющей такое построение, не вызывает осложнений.
Рядом Фурье периодической с периодом 2l функции f(x) называется тригонометрический ряд вида











1 0
sin cos
2 1
n
n
n
l
x
n
b
l
x
n
a
a


в котором коэффициенты вычисляются по следующим формулам



l
l
dx
x
f
l
a
)
(
2 0
,



l
l
n
dx
l
x
n
x
f
l
a

cos
)
(
1
,



l
l
n
dx
l
x
n
x
f
l
b

sin
)
(
1
Создадим процедуру разложения в ряд Фурье алгебраического выражения.
Ее параметрами являются само выражение, имя независимой переменной, по которой выражение раскладывается в ряд, значение полупериода l и количество удерживаемых членов.
FourierSeries:=proc(f::algebraic,x::name, l::anything,n::nonnegint)
(int(f,x=-l..l)+ sum(int(f*cos(k*Pi*x/l),x=-l..l)*cos(k*Pi*x/l),k=1..n)+ sum(int(f*sin(k*Pi*x/l),x=-l..l)*sin(k*Pi*x/l),k=1..n))/l; end proc:

2
Пример 1. Разложим на отрезке
]
,
[



в ряд Фурье функцию f(x)=x.
> ser1:=normal(FourierSeries(f,x,Pi,6));
> plot(ser1,x=-3*Pi..3*Pi,thickness=2,color=BLACK);
Поскольку функция нечетная, то в разложении нет косинусов. По графику ясно, что удерживаемых в ряде членов недостаточно для приемлемого приближения пилообразной функции (на отрезке
]
,
[



функция равна x и вне него периодически повторяет свои значения). Увеличение числа удерживаемых членов ряда позволяет повысить точность приближения.
Отметим, однако, что в окрестности точек разрыва ряд сильно осциллирует и его значения сильно отличаются от значений исходной функции.
Пример 2. Разложим с отрезка [-2,2] в ряд Фурье пилообразную функцию следующего вида f:=-PR(x,-2,1)+2*PR(x,-1,2)-PR(x,1,1): plot(f,x=-3..3,thickness=2,color=BLACK);
Ее уравнение и ряд имеют следующий вид f; ser2:=FourierSeries(f,x,2,6);
Сама функция довольно хорошо приближается рядом Фурье. На следующем рисунке слева показан график ряда при 6 удерживаемых членах (предыдущая формула). На рисунке справа показан график ряда при 30-ти удерживаемых членах.
> plot(ser2,x=-4..4,thickness=2,color=BLACK);





2
( )
sin x
(
)
sin 2 x
2 3
(
)
sin 3 x
1 2
(
)
sin 4 x
2 5
(
)
sin 5 x
1 3
(
)
sin 6 x




1 2

x
2

x
1

x
1 1
2

x
2
:=
ser2


8






sin

x
2

2 8
9






sin
3

x
2

2 8
25






sin
5

x
2

2

3
Увеличения числа членов ряда делает приближение еще лучшим.
Если функцию, разлагаемую в ряд Фурье, надо определять в виде процедуры, то в нашу процедуру разложения в ряд надо внести поправки. В этом случае она принимает вид
FourierSeries1:=proc(f::procedure,x::name, l::anything,n::nonnegint)
(int(f(x),x=-l..l)+ sum(int(f(x)*cos(k*Pi*x/l),x=-l..l)*cos(k*Pi*x/l),k=1..n)+ sum(int(f(x)*sin(k*Pi*x/l),x=- l..l)*sin(k*Pi*x/l),k=1..n))/l; end proc:
Здесь первый параметр должен быть типа процедуры и в теле вместо выражения f следует использовать f(x).
Пример 3. Разложим в ряд Фурье с отрезка
]
,
[



функцию
x
x
f

)
(
f:=proc(x) abs(x) end proc:
ser3:=collect(FourierSeries1(f,x,Pi,10),cos); plot(ser3,x=-2*Pi..2*Pi,thickness=2,color=BLACK);
Как видим, уже при 10-ти удерживаемых членах график ряда хорошо аппроксимирует на отрезке
]
,
[



функцию
x
x
f

)
(
. Отметим также, что в разложении в ряд четной функции
x
x
f

)
(
присутствуют только косинусы.
Пример 4. Разложим функцию
2 1
)
(
x
x
f


в ряд Фурье на отрезке [-1,1]. f:=x->1-x^2: ser5:=collect(FourierSeries1(f,x,1,6),cos); plot(ser5,x=-3..3,thickness=2,color=BLACK,scaling=CONSTRAINED);
Разложение четной функции и здесь дает ряд только по косинусам.






4 81
(
)
cos 9 x

4 49
(
)
cos 7 x

4
( )
cos x

4 9
(
)
cos 3 x

4 25
(
)
cos 5 x



4
2. Решение задачи Коши разложением в степенной ряд.
Применение рядов для решения ДУ рассмотрим на примере нелинейного уравнения второго порядка вида
)
,
,
(
y
y
x
f
y



удовлетворяющего начальным условиям
0 0
0 0
)
(
,
)
(
y
x
y
y
x
y




Предположим, что решение поставленной задачи существует и будем его искать в форме ряда Тейлора функции y(x)









2 0
0 0
0 0
)
)(
(
2 1
)
)(
(
)
(
)
(
x
x
x
y
x
x
x
y
x
y
x
y
Для построения решения нужно знать производные в точке x
0
порядка 1, 2, 3,
…, но это можно сделать с помощью самого уравнения и его начальных условий. Действительно, из начальных условий мы сразу имеем значение функции и ее первой производной в точке x
0
. Подставляя их в уравнение, сразу находим значение
2-й производной
)
,
,
(
)
(
0 0
0 0
y
y
x
f
x
y



Дифференцируя обе части исходного уравнения по аргументу x, получим выражение для 3-й производной
y
y
y
y
x
f
y
y
y
y
x
f
x
y
y
x
f
y
















)
,
,
(
)
,
,
(
)
,
,
(
Подставляя в правую часть x=x
0
и значения предыдущих производных в этой точке, находим значение 3-й производной
)
(
0
x
y

в точке x=x
0
Дифференцируя последнее выражение еще раз, и выполняя подстановку известных производных, находим значение 4-й производной в точке x=x
0
Для определения производных функции y(x) в точке x
0
до требуемого порядка n дифференцирование и подстановку выполняем необходимое число раз. Затем найденные значения производных подставляем в представление решения в форме ряда Тейлора. Для тех значений x для которых этот ряд сходится, он будет представлять искомое решение задачи Коши.
Для решения задачи Коши
)
,
,
(
y
y
x
f
y



,
0 0
0 0
)
(
,
)
(
y
x
y
y
x
y




используем следующую процедуру (см. [1], стр. 450).
DiffSer:=proc(f::anything,x::name,y::name,y1::name,\ x0::numeric,y0::numeric,y01::numeric,n::integer) local p,der,i,j,s: m0||1:=y01: m0||2:=subs(x=x0,y=y0,y1=y01,f); der:=subs(y1=m1,f); p:=y0+m0||1*(x-x0)+m0||2*(x-x0)^2/2: for i from 3 to n do s:=0: for j from 1 to i-2 do s:=s+diff(der,m||j)*m||(j+1); end do: der:=diff(der,x)+diff(der,y)*m||1+s; m0||i:=eval(der,{x=x0,y=y0,seq(m||k=m0||k,k=1..i-1)}); p:=p+(m0||i)*(x-x0)^i/i!; end do; end proc:

5
Аргументами процедуры являются выражение для правой части ДУ (f), имена независимой переменной (x), неизвестной функции (y) и ее производной (y1), а также параметры x0, y0, y01, определяющие условия задачи Коши. Порядок усечения ряда n, задаваемый последним параметром, указывает, что последний член ряда будет иметь вид
n
n
x
x
a
)
(
0

Операторы расположенные до первого цикла for формируют первые 3 члена ряда. В цикле for вычисляются оставшиеся производные и их значения в точке x=x
0
, а также формируется ряд. В процедуре для всех появляющихся при последовательном дифференцировании производных решения используются обозначения mN, где N – порядок производной, а для их значений в точке x
0
используются переменные m0N.
Пример 1. Решить следующую задачу Коши:
,
, s:=DiffSer(-y*x^2,x,y,y1,0,1,0,20);
Выбранная задача имеет точное решение. dsolve({diff(y(x),x$2)=-y(x)*x^2,y(0)=1,D(y)(0)=0},y(x)); s1:=rhs(%);
Для сравнения построим графики точного и приближенного решений. На следующем графике слева показано точное решение и приближенное до 20 –
й степени, а справа до 100-й степени..
plot([s,s1],x=0..5,-
1..1,color=BLACK,thickness=2,linestyle=[1,4]);
Как видим, увеличение числа членов ряда расширяет интервал аппроксимации решения.
Пример 2. Исследуем колебания плоского математического маятника – точка массы m подвешена на конце нити длины l и находится в однородном поле тяжести. Колебания маятника описываются уравнением
0
sin
2





mgl
ml
, где
)
(t

- угол отклонения маятника от положения равновесия. Перепишем это уравнение в виде


sin
l
g
tt



. Это вид ДУ, который можно решать с помощью процедуры DiffSer.



2
x
2
y

y x
2

( )
y 0 1

d
d
x
( )
y 0 0
:=
s





1 1
12
x
4 1
672
x
8 1
88704
x
12 1
21288960
x
16 1
8089804800
x
20

( )
y x

1 2







3 4
x






BesselJ
,
1 4
x
2 2
1 2







3 4
x






BesselY
,
1 4
x
2 2

6
Для примера выберем m=1, g=1, а начальные условия возьмем следующего вида:
8 1
)
0
(
,
0
)
0
(





т.е. задана начальная скорость и в нулевой момент времени маятник находится в положении равновесия.
Решим эту задачу Коши с помощью процедуры
DiffSer
> s:=DiffSer(-sin(y),x,y,y1,0,0,1/8,14);
> plot(s,x=0..Pi+1,color=BLACK,thickness=2, linestyle=[1,4],numpoints=1000);
Если начальная скорость не слишком велика, то колебания должны быть периодическими. Мы видим, что наше решение через некоторый промежуток времени принимает нулевое значение, а значит, маятник возвращается в положение равновесия. Вероятно о дальнейшем поведении решения наш ряд ничего не говорит, поскольку начинает сказываться погрешность приближения.
Пример 3. Решим то же ДУ с другими начальными условиями
0
)
0
(
,
8 1
)
0
(





, т.е. задана начальное отклонение относительно положения равновесия, а начальная скорость равна 0.
> s:=DiffSer(-sin(y),x,y,y1,0,1/8,0,9): evalf[5](s);
> plot(s,x=0..Pi+2,color=BLACK,thickness=2, linestyle=[1,4],numpoints=1000);
Из графика видно, что материальная точка возвращается в положение равновесия, затем удаляется от него в противоположном направлении. Для больших моментов времени полученное решение применять нельзя.
Построить решение в форме ряда можно и другим способом. В этом способе решение y(x) представляется в виде степенного ряда с неопределенными коэффициентами









3 0
3 2
0 2
0 1
0
)
(
)
(
)
(
)
(
x
x
a
x
x
a
x
x
a
a
x
y
s
1 8
x
1 48
x
3 13 12288
x
5 4801 165150720
x
7 19531 21743271936
x
9




:=
38052343 765363172147200
x
11 4421087563 1528277182143528960
x
13






0.12500 0.062335 x
2 0.0051542 x
4 0.00016239 x
6 0.14468 10
-5
x
8

7
Этот ряд подставляется в левую и правую части уравнения
)
,
,
(
y
y
x
f
y



, причем все функции правой части также раскладываются в степенные ряды.
Далее выполняются необходимые вычисления с рядами правой части для представления ее в виде одного степенного ряда. Приравнивая в левой и правой части коэффициенты при одинаковых степенях (x-x
0
)
n
, получим систему линейных уравнений относительно неизвестных коэффициентов ряда a
n
В работе [1] стр.453 приведена процедура решения задачи Коши для
ДУ 2-го порядка, имеющей вид
)
,
,
(
y
y
x
f
y



,
0 0
0 0
)
(
,
)
(
y
x
y
y
x
y




, реализующая описанный алгоритм. Вот ее код. intDiff:=proc(f::anything,x::name,y::name,y1::name,\ x0::numeric,y0::numeric,y01::numeric,n::integer) local p,pL,i,s,a,result: p:=y0+y01*(x-x0): result:=p+sum(a[i]*(x-x0)^i,i=2..n); p:=p+sum(a[i]*(x-x0)^i,i=2..n+2); pL:=series(subs(y=p,y1=diff(p,x),f),x=x0,n+1); for i from 0 to n do eq||i:=coeff(diff(p,x$2)-convert(pL,polynom),x,i)=0; end do: s:=solve({seq(eq||i,i=0..n)},{seq(a[i],i=2..n+2)}); assign(s); eval(result); end proc:
Аргументы этой процедуры имеют тот же смысл, что и аргументы предыдущей процедуры. В теле процедуры введен промежуточный ряд, имеющий слагаемых на 2 больше, чем требуемое число членов, т.к. его приходится дважды дифференцировать. Отметим, что эта процедура работает быстрее предыдущей.
Пример 4. Решим задачу Коши
0





xy
y
y
x
,
0
)
0
(
,
1
)
0
(



y
y
. Для решения это уравнение с помощью разработанной процедуры, преобразуем его к виду
y
x
y
y





. Отметим, что функция правой части при x=0 не определена. Поэтому применение процедуры реализующей первый способ построения решения в виде ряда не представляется возможным, т.к. ее алгоритм предполагает вычисление правой части в точке x=0. Но в методе поиска неопределенных коэффициентов, реализованном в последней процедуре, значение правой части в этой точке не вычисляется.
> s:=intDiff(-y1/x-y,x,y,y1,0,1,0,12);
> s1:=dsolve({x*diff(y(x),x$2)+diff(y(x),x)+x*y(x)=0, y(0)=1,D(y)(0)=0},y(x));
> plot([s,rhs(s1)],x=0..7,color=BLACK,thickness=2, linestyle=[1,4],numpoints=1000);






1 1
4
x
2 1
64
x
4 1
2304
x
6 1
147456
x
8 1
14745600
x
10 1
2123366400
x
12
:=
s1

( )
y x
(
)
BesselJ ,
0 x

8
На графике приведено точное решение задачи и приближенное в виде ряда.
Расхождение решения в форме ряда с точным решением может происходить по двум причинам. Первая состоит в том, что мы не можем выбрать достаточное количество членов ряда для получения удовлетворительной точности на широком интервале изменения аргумента.
Это ограничение связано с ограничением по времени выполнения программы и по объему используемой ею памяти. Вторая причина состоит в том, что большинство рядов имеют конечный радиус сходимости и использование ряда для приближения решения на более широком интервале принципиально невозможно. Выход из создавшегося положения состоит в разбиении интервала поиска решения на отрезки такой длины, в пределах которых степенные ряды будут сходиться и давать приемлемое приближение решения. Вначале строится решение s
1
(x) в форме ряда Тейлора на отрезке
[x
0
,x
1
], где мы знаем, что ряд хорошо аппроксимирует решение и сходится.
Затем находим значение этого ряда и его производной в точке x
1
и используем эти значения в качестве начальных значений для построения нового ряда s
2
(x) разложения решения в окрестности точки x
1
по степеням (x-
x
1
)
n
. Построенный ряд используем для определения начальных условий Коши в точке x=x
2
и определяем третий ряд аппроксимирующий решение на следующем отрезка поиска решения и т.д.
Для иллюстрации этого подхода рассмотрим ДУ описывающее колебания математического маятника.
Пример 5. Для нашей задачи возьмем начальные условия вида
0
)
0
(
,
8 1
)
0
(



y
y
и для построения решения в форме ряда используем процедуру intDiff. s1:=intDiff(-sin(y),x,'y','y1',0,0,1/8,20): evalf[3](s1); plot(s1,x=0..Pi+1,color=BLACK,thickness=2, linestyle=[1,4],numpoints=1000);
0.125 x
0.0208 x
3 0.00106 x
5 0.0000291 x
7 0.898 10
-6
x
9




0.497 10
-7
x
11 0.289 10
-8
x
13 0.146 10
-9
x
15 0.697 10
-11
x
17




0.351 10
-12
x
19


9
Полученный ряд дает удовлетворительное решение на отрезке [0,3].
Используем x=3 в качестве начальной точки для решения новой задачи
Коши. Подставим x=3 в найденный ряд и его производную по x для получения новых значений
0 0
)
0
(
,
)
0
(
y
y
y
y




. Опять решаем задачу Коши используя процедуру intDiff. x0:=3: y0:=subs(x=x0,s1): y1:=subs(x=x0,diff(s1,x$1)): s2:=intDiff(-sin(y),x,'y','y1',x0,y0,y1,14): evalf[3](s2);
Из двух рядов создаем кусочно-непрерывную функцию и строим ее график.
Sser:=x->piecewise(x<=x0,s1,x>x0,s2): plot(Sser(x),x=0..4,color=BLACK,thickness=2, linestyle=[1,4],numpoints=1000);
Полученная функция уже хорошо приближает решение на отрезке [0,6], где она похожа на синусоиду. Мы это знаем т.к. решение должно быть нечетным и периодическим (задача о колебании маятника). Интересно отметить, что построенная функция непрерывна вместе со своей 1-й производной.
Действительно, 1-я производная первого ряда в точке x=3 использована в качестве первой производной второго ряда как начальное значение. Вторая производная в точке стыковки для обеих функций, как правило, будет различной. Для первого ряда она вычисляется подстановкой x=x
1
во 2-ю производную первого ряда, а для второй вычисляется по формуле
)
,
,
(
)
(
1 1
1 1
y
y
x
f
x
y



Для реализации метода решения задачи Коши ДУ 2-го порядка разложением в степенной ряд с разбиением на отрезки мы создали специальную процедуру.
StepSeries:=proc(f::anything,x::name,y::name,y1::name,\ x0::numeric, y0::numeric, y01::numeric,\ n::integer,xend::numeric, dlt::numeric) local i,nn,DiffInner;
# Внутренняя процедура разложения решения в ряд
# в окрестности 0
DiffInner:=proc(f::anything,x::name,y::name,y1::name,\ y0::Or(numeric,realcons),\ y01::Or(numeric,realcons),n::integer)
0.392 0.125 x
0.00895 (
)

x
3.
2 0.0208 (
)

x
3.
3 0.000758 (
)

x
3.
4




0.00106 (
)

x
3.
5 0.0000291 (
)

x
3.
6 0.0000287 (
)

x
3.
7



0.115 10
-5
(
)

x
3.
8 0.857 10
-6
(
)

x
3.
9 0.771 10
-7
(
)

x
3.
10



0.461 10
-7
(
)

x
3.
11 0.530 10
-8
(
)

x
3.
12 0.265 10
-8
(
)

x
3.
13



0.306 10
-9
(
)

x
3.
14


10 local p,pL,i,s,a,result: p:=y0+y01*x: result:=p+sum(a[i]*x^i,i=2..n); p:=p+sum(a[i]*x^i,i=2..n+2); pL:=series(subs(y=p,y1=diff(p,x),f),x=0,n+1); for i from 0 to n do eq||i:=coeff(diff(p,x$2)-convert(pL,polynom),x,i)=0; end do: s:=solve({seq(eq||i,i=0..n)},{seq(a[i],i=2..n+2)}); assign(s); eval(result); end proc: nn:=floor((xend-x0)/dlt): # число точек стыковки
# Второе число X0i,Y0i,Y1i указывает номер шага i
X0||0:=x0: Y0||0:=y0: Y1||0:=y01: s||1:=DiffInner(f,x,y,y1,Y00,Y10,n): for i from 1 to nn-1 do
X0||i:=x0+dlt*i:
Y0||i:=evalf(subs(x=dlt,s||i));
Y1||i:=evalf(subs(x=dlt,diff(s||i,x$1))); s||(i+1):=DiffInner(f,x,y,y1,Y0||i,Y1||i,n): end do: piecewise(seq(op([x<=X0||i,subs(x=x-X0||(i-
1),s||i)]),i=1..nn-1),\ subs(x=x-X0||(nn-1),s||nn)); end proc:
Аргументами процедуры являются выражение для правой части ДУ (f), имена независимой переменной (x), неизвестной функции (y) и ее производной (y1), а также параметры x0, y0, y01, определяющие условия задачи Коши. Порядок усечения рядов задаваемый параметром n, указывает, что их последние члены будут иметь вид
n
n
x
x
a
)
(
0

. Эти аргументы в точности те же, что и у процедуры DiffSer. Последние два параметра определяют конец интервала поиска решения xend (оно ищется на интервале
[x0, xend]) и длину отрезков на которые разбивается этот интервал – dlt. Это означает, что будут построено






dlt
x
xend
N
0
степенных рядов, которые будут стыковаться в точках
N
k
dlt
k
x
x
k
,
,
2
,
1
,
0





. Внутренняя процедура DiffInner
идентична процедуре DiffSer и всегда вызывается со значением x0=0, поэтому DiffInner не имеет этого аргумента. Результатом работы процедуры является piecewise функция, составленная из построенных рядов для которых выполняется преобразование к общей системе координат
(координатная система первого отрезка).
Работу процедуры проиллюстрируем на примере уравнения колебания математического маятника. Мы не будем приводить формулы полученной
piecewise функции из-за их длины.