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

11
Пример 6. Для маятника задано начальное отклонение, а начальная скорость равна 0. st:=time(): x0:=0: y0:=1/8: y1:=0: r1:=StepSeries(-sin(y),x,'y','y1',x0,y0,y1,8,15,3): time()-st; plot(r1,x=-5..11.5,color=BLACK,thickness=2,numpoints=500);
На следующем рисунке слева приведен график решения, составленного из многочленов до 8-й степени, а справа –14-й. Оба решения составлены из 5-ти рядов построенных на отрезках [0,3], [3,6], [6,9], [9,12], [12,15]. На левом рисунке заметно некоторое отличие 3-й волны графика от 2-х предыдущих, поэтому для улучшения решения пришлось повысить степень многочленов до 14-ти.
Выбрав другое начальное условие (начальное отклонение 0, начальная скорость не 0) получим следующее решение x0:=0: y0:=0: y1:=1/8: r1:=StepSeries(-sin(y),x,'y','y1',x0,y0,y1,20,15,3): plot(r1,x=-4.1..16.5,color=BLACK,thickness=2,numpoints=500);
Здесь для построения приемлемого решения пришлось строить многочлены
20-й степени. Полученное решение можно применять на интервале несколько более широком, чем отрезок заданный в аргументе процедуры, поскольку первый ряд применим левее начальной точки x0, а правый может давать приближение на отрезке несколько более широком, чем dlt. Если мы правильно выбрали параметр dlt, то решение в виде составного ряда будет всегда приемлемым на отрезке [x0-dlt, xend].
Если начальная скорость достаточно высока (в нашем примере больше
2), то движение маятника перестает быть периодическим (маятник доходит до верхней точки и переваливает через нее). В этом случае угол отклонения маятника возрастает монотонно. st:=time(): x0:=0: y0:=0: y1:=2.05: r1:=StepSeries(-sin(y),x,'y','y1',x0,y0,y1,12,9,1): time()-st; plot(r1,x=-1..10,color=BLACK,thickness=2,numpoints=500);

12
На следующем рисунке слева показан график решения при нулевом начальном смещении и начальной скорости y1=2.05, а на правом – при начальной скорости равной y1=3. Если начальная скорость достаточно велика, то кривизны графика решения почти незаметно. Легко видеть, что вне отрезка [0,9] построенные решения неприменимы.
Отметим, что для построения последних решений длину отрезков, на которых построены ряды, пришлось брать равной единице, т.е. наши решения составлены из 9 многочленов 12-й степени. Если бы мы не использовали процедуру сшивания рядов, то не смогли бы построить хорошее решение на таком широком интервале изменения аргумента.
Отметим также, что периодическое решение, соответствующее колебательному процессу достаточно построить на отрезке равном периоду колебаний, а на более широкий интервал времени решение мы можем продолжить с помощью нашей методики периодического продолжения функций.
3.
Метод Галеркина решения краевых задач ОДУ
Общая краевая задача для ОДУ 2-го порядка формулируется следующим образом - требуется найти решение ДУ вида
b
x
a
x
f
u
x
q
u
x
p
u
Lu








),
(
)
(
)
(
, удовлетворяющее краевым условиям
)
(
)
(
,
)
(
)
(
1 1
1 1
0 0
0 0














a
u
a
u
u
l
a
u
a
u
u
l
где p, q и f являються непрерывными функциями переменой x на интервале
[a,b], а

i
,

i
,

i
– некоторые числа.
Для нахождения приближенного решения краевой задачи на отрезке
[a,b] задают некоторую линейно независимую систему дважды непрерывно дифференцируемых функций
,...
,...,
,
1 0
n



таких, что функция
0

удовлетворяет заданным краевым условиям
1 0
1 0
0 0
,






l
l
, а остальные функции системы
i

при i=1,2,… удовлетворяют однородным краевым условиям
0
,
0 1
0


i
i
l
l


. Эта система функций называется базисной.
Приближенное решение краевой задачи ищется в виде
)
(
)
(
)
(
)
(
1 1
0
x
a
x
a
x
x
y
n
n
n








13 с неизвестными коэффициентами a
i
, i=1,2,…,n., которые определяются из выполнения некоторого критерия. Отметим, что в силу выбора базисных функций функция y
n
(x) будет точно удовлетворять краевым условиям.
В методе Галеркина в качестве критерия для определения неопределенных коэффициентов выбирается требование ортогональности базисных функций
n



,..,
,
2 1
к функции невязки
)
,...,
,
(
1
n
a
a
x

являющейся разностью левой и правой части ОДУ при подстановке в него приближенного решения y
n
(x).








n
k
k
k
n
n
x
L
a
x
f
x
L
x
f
Ly
a
a
x
1 0
1
)
(
)
(
)
(
)
(
)
,...,
,
(


Под ортогональностью понимается равенство нулю скалярного произведения невязки и базисных функций
0
)
(
)
,...,
,
(
1



b
a
i
n
dx
x
a
a
x

. Это требование приводит к линейной системе алгебраических уравнений для определения неизвестных коэффициентов a
i
приближенного решения y
n
(x).
)
,
(
)
,
(
)
,
(
)
,
(
)
,
(
)
,
(
)
,
(
)
,
(
)
,
(
0 1
1 2
0 2
2 1
1 1
0 1
1 1
1
n
n
n
n
n
n
n
n
n
L
f
L
a
L
a
L
f
L
a
L
a
L
f
L
a
L
a






























Решение этой системы всегда существует и единственно в силу линейной независимости базисных функций.
Вот процедура построения приближенного решения краевой задачи
ОДУ второго порядка, реализующая метод Галеркина (см. [1] стр.474).
Galerkin:=proc(L::procedure,f::algebraic,phi::procedure, x::name,a::numeric,b::numeric,n::numeric) local y,i,eq,A: y:=phi(0,x)+sum(A[i]*phi(i,x),i=1..n): for i to n do eq||i:=sum(A[k]*int(L(phi(k,x),x)*phi(i,x),x=a..b),k=1..n)=\ int((f-L(phi(0,x),x))*phi(i,x),x=a..b); end do: eq:=solve({seq(eq||i,i=1..n)},{seq(A[i],i=1..n)}); assign(eq); simplify(eval(y)); end proc:
Параметрами этой процедуры является дифференциальный оператор L, правая часть уравнения f, процедура вычисления базисных функций phi.
Параметр x определяет имя независимой переменной, a,b – координаты граничных точек, n – число удерживаемых членов в приближенном решении. Дифференциальный оператор L передается в виде процедуры с двумя параметрами – алгебраическим выражением и именем независимой переменной, по которой выполняется дифференцирование. Параметр- процедура phi вычисляет n-ю базовую функцию
n

, и ее параметрами являются номер базовой функции и имя независимой переменной.

14
Пример1. Решим методом Галеркина краевую задачу
0
)
1
(
,
0
)
0
(
,
1 0
,










u
u
x
x
u
u
Lu
Выбираем базисные функции
,...
2
,
1
),
1
(
)
(
,
0
)
(
0




i
x
x
x
x
i
i


. Функция
)
(
0
x

удовлетворяет граничным условиям краевой задачи, а функции
)
(x
i

однородным граничным условиям.
Создадим необходимые процедуры вычисления дифференциального оператора, базовых функций и выражения правой части.
> L1:=proc(y,x) diff(y,x$2)+y; end proc:
> f:=-x:
> phi:=proc(n,x) if n=0 then 0 else x^n*(1-x) end if; end proc:
Построим приближенное решение
> res1:=Galerkin(L1,f,phi,x,0,1,1); res2:=Galerkin(L1,f,phi,x,0,1,2); res3:=Galerkin(L1,f,phi,x,0,1,3);
Найдем точное решение этой краевой задачи
> exact:=dsolve({diff(y(x),x$2)+y(x)=-x,y(0)=0,y(1)=0},y(x));
Для сравнения построим графики полученных решений
> plot([res1,res2,res3,rhs(exact)],x=0..1,color=black, thickness=2,linestyle=[4,7,1]);
Видно, что приближенные решения, в которых указаны две базовые функции
(и три) на графике сливаются с точным.
Пример 2. Решим то же ДУ с другими краевыми условиями
0
)
1
(
,
1
)
0
(


u
u
Для этого надо выбрать другую функцию
x
x


1
)
(
0

и, следовательно, изменить процедуру phi.
> phi:=proc(n,x) x^n*(1-x);
:=
res1

5 x (
)
 
1
x
18
:=
res2


71 369
x
8 369
x
2 7
41
x
3
:=
res3



13811 73554
x
469 73554
x
2 2667 12259
x
3 7
299
x
4
:=
exact

( )
y x

( )
sin x
( )
sin 1
x

15 end proc:
> res1:=Galerkin(L1,f,phi,x,0,1,1); res2:=Galerkin(L1,f,phi,x,0,1,2); res3:=Galerkin(L1,f,phi,x,0,1,3); exact:=dsolve({diff(y(x),x$2)+y(x)=-x,y(0)=1,y(1)=0},y(x));
Построим графики полученных решений
График уже первого приближения сливается с графиком точного решения.
Пример 3. Метод Галеркина имеет смысл применять для задач, не имеющих точного решения. Пусть требуется решить следующую краевую задачу
0
)
1
(
,
0
)
1
(
,
1 1
,
1
)
1
(
2













u
u
x
u
x
u
Lu
Выбираем базисные функции
,...
2
,
1
),
1
(
)
(
,
0
)
(
2 2
2 0





i
x
x
x
x
i
i


. Функция
)
(
0
x

удовлетворяет граничным условиям краевой задачи, а функции
)
(x
i

однородным граничным условиям.
Создадим необходимые процедуры вычисления дифференциального оператора, базовых функций и выражения правой части
> L1:=proc(y,x) diff(y,x$2)+(1+x^2)*y; end proc:
> f:=-1:
> phi:=proc(n,x) if n=0 then 0 else x^(2*n-2)*(1-x^2) end if: end proc:
Решение
Находим приближенные решения при n=1, 2, 3.
> res1:=Galerkin(L1,f,phi,x,-1,1,1); res2:=Galerkin(L1,f,phi,x,-1,1,2); res3:=Galerkin(L1,f,phi,x,-1,1,3);
> plot([res1,res2,res3],x=-1..1,color=black, thickness=2,linestyle=[4,7,1]);
:=
res1
 
1 4
9
x
5 9
x
2
:=
res2
 
1 4
9
x
5 9
x
2
:=
res3




1 407 897
x
448 897
x
2 28 299
x
3 14 299
x
4
:=
exact

( )
y x



( )
sin x (
)

( )
cos 1 1
( )
sin 1
( )
cos x
x
:=
res1

35 38 35 x
2 38

16
Поскольку точного решения не существует, то для оценки погрешности приближенных решений мы построим график невязок.
> plot([L1(res1,x)-f,L1(res2,x)-f,L1(res3,x)-f],x=-
1..1,color=black,thickness=2,linestyle=[4,7,1]);
Невязка третьего приближения отображена сплошной линией и практически полностью сливается с осью x.Для оценки точности приведем ее график отдельно.
Максимальная невязка для третьего приближения составляет 0.014
Метод Галеркина весьма эффективен для решения краевых задач ОДУ.
Однако при его применении возникает трудность выбора линейно независимой системы базисных функций, удовлетворяющей однородным краевым условиям. В связи с этим отметим, что получивший в последнее время широкое распространение метод конечных элементов является дальнейшим развитием метода Галеркина. Фактически МКЭ имеет специфические особенности выбора базисной системы функций, а в остальном совпадает с методом Галеркина.
Литература.
1. А.. Матросов. “Maple 6. Решение задач высшей математики и механики”,
БХВ, Петербург, 2001г.
:=
res2


3969 4252 1050 1063
x
2 231 4252
x
4
:=
res3



5148231 5523436 62073 64226
x
2 26169 5523436
x
4 81939 2761718
x
6