Задача Штурма-Лиувилля

Дисциплина:	Раздел:	Выполнил:
Методы конечных элементов	Одномерные конечные элементы. Уравнение Штурма-Лиувилля	ЧГУ, матфак, 4 курс Осипов.П.Г

www.domath.ru

`>`	#Подключим пакет графики plots.

`>`	restart: with(plots):

`>`	# Начальные данные, [a,b] - отрезок, d,e - граничные условия, f,p и q исходные функции.

`>`	a:=3:b:=7: d:=0.49:e:=-0.12: f:=x -> ln(x): p:=x -> arctan(x): q:=x-> ln(ln(x)):

`>`	#Случай для 10 узлов

> n:=10:t[0]:=a:c[0]:=d:

#В даннном цикле зададим равномерную сетку.
for i from 1 to n
do
t[i]:=t[i-1]+(b-a)/n
end do:
#Зададим кусочную функцию по методу.
func_piece(x,0):= piecewise( (x>=t[0]) and (x<t[1]),
(t[1]-x)/(t[1]-t[0]),
(x>=t[1]),
0):
func_piece(x,n):= piecewise( (x>=t[n-1]) and (x<=t[n]),
(x-t[n-1])/(t[n]-t[n-1]),
(x<t[n-1]),
0):

for i from 1 to n-1
do
func_piece(x,i):= piecewise( (x>=t[i-1]) and (x<t[i]),

(x-t[i-1])/(t[i]-t[i-1]),
(x>=t[i]) and (x<t[i+1]),
(t[i+1]-x)/(t[i+1]-t[i]),
(x<t[i-1]) and (x>=t[i+1]), 0)
end do:

#Зададим систему алгебраических уравнений относительно c[i}.
r[1]:=evalf(

c[1]*(Int(p(x)*diff(func_piece(x,1),x)^2,x=t[0]..t[2])+
Int(q(x)*func_piece(x,1)^2,x=t[0]..t[2]))+

c[2]*(Int(p(x)*diff(func_piece(x,2),x)*diff(func_piece(x,1),x),x=t[1]..t[2])+
Int(q(x)*func_piece(x,2)*func_piece(x,1),x=t[1]..t[2]))-

Int(func_piece(x,1)*f(x),x=t[0]..t[2])+ d*(Int(p(x)*diff(func_piece(x,0),x)*diff(func_piece(x,1),x),x=t[0]..t[1])+

Int(q(x)*func_piece(x,0)*func_piece(x,1),x=t[0]..t[1]))):

r[n]:=evalf(

c[n-2]*(Int(p(x)*diff(func_piece(x,n-2),x)*diff(func_piece(x,n-1),x),x=t[n-2]..t[n-1])+
Int(q(x)*func_piece(x,n-2)*func_piece(x,n-1),x=t[n-2]..t[n-1]))+

c[n-1]*(Int(p(x)*diff(func_piece(x,n-1),x)^2,x=t[n-2]..t[n])+
Int(q(x)*func_piece(x,n-1)^2,x=t[n-2]..t[n]))-
Int(func_piece(x,n-1)*f(x),x=t[n-2]..t[n])+

e*(Int(p(x)*diff(func_piece(x,n-1),x)*diff(func_piece(x,n),x),x=t[n-1]..t[n])
+Int(q(x)*func_piece(x,n-1)*func_piece(x,n),x=t[n-1]..t[n]))):

for k from 2 by 1 to n-1
do

r[k]:=evalf(

c[k-1]*(Int(p(x)*diff(func_piece(x,k-1),x)*diff(func_piece(x,k),x),x=t[k-1]..t[k]) +Int(q(x)*func_piece(x,k-1)*func_piece(x,k),x=t[k-1]..t[k]))+

c[k]*(Int(p(x)*diff(func_piece(x,k),x)^2,x=t[k-1]..t[k+1])
+Int(q(x)*func_piece(x,k)^2,x=t[k-1]..t[k+1]))+

c[k+1]*(Int(p(x)*diff(func_piece(x,k+1),x)*diff(func_piece(x,k),x),x=t[k]..t[k+1])
+Int(q(x)*func_piece(x,k)*func_piece(x,k+1),x=t[k]..t[k+1]))-
Int(func_piece(x,k)*f(x),x=t[k-1]..t[k+1])):

end do:

for i from 1 to n
do
c[i]:= solve(r[i], c[i])
end do:

y_10(x):=c[0]*func_piece(x,0)+sum(c[j]*func_piece(x,j),j=1..n):

`>`	#Случай для 40 узлов.

> m:=40:t1[0]:=a:c1[0]:=d:
#Построим сетку.
for i from 1 to m
do
t1[i]:=t1[i-1]+(b-a)/m
end do:
#Зададим кусочные функции согласно методу.
func_piece1(x,0):= piecewise( (x>=t1[0]) and (x<t1[1]),
(t1[1]-x)/(t1[1]-t1[0]),
(x>=t1[1]),
0):
func_piece1(x,m):= piecewise( (x>=t1[m-1]) and (x<=t1[m]),
(x-t1[m-1])/(t1[m]-t1[m-1]),
(x<t1[m-1]),
0):

for i from 1 to m-1
do
func_piece1(x,i):= piecewise( (x>=t1[i-1]) and (x<t1[i]),

(x-t1[i-1])/(t1[i]-t1[i-1]),
(x>=t1[i]) and (x<t1[i+1]),
(t1[i+1]-x)/(t1[i+1]-t1[i]),
(x<t1[i-1]) and (x>=t1[i+1]), 0)
end do:

r1[1]:=evalf(

c1[1]*(Int(p(x)*diff(func_piece1(x,1),x)^2,x=t1[0]..t1[2])+
Int(q(x)*func_piece1(x,1)^2,x=t1[0]..t1[2]))+

c1[2]*(Int(p(x)*diff(func_piece1(x,2),x)*diff(func_piece1(x,1),x),x=t1[1]..t1[2])+
Int(q(x)*func_piece1(x,2)*func_piece1(x,1),x=t1[1]..t1[2]))-

Int(func_piece1(x,1)*f(x),x=t1[0]..t1[2])+ d*(Int(p(x)*diff(func_piece1(x,0),x)*diff(func_piece1(x,1),x),x=t1[0]..t1[1])+

Int(q(x)*func_piece1(x,0)*func_piece1(x,1),x=t1[0]..t1[1]))):

r1[m]:=evalf(

c1[m-2]*(Int(p(x)*diff(func_piece1(x,m-2),x)*diff(func_piece1(x,m-1),x),x=t1[m-2]..t1[m-1])+
Int(q(x)*func_piece1(x,m-2)*func_piece1(x,m-1),x=t1[m-2]..t1[m-1]))+

c1[m-1]*(Int(p(x)*diff(func_piece1(x,m-1),x)^2,x=t1[m-2]..t1[m])+
Int(q(x)*func_piece1(x,m-1)^2,x=t1[m-2]..t1[m]))-
Int(func_piece1(x,m-1)*f(x),x=t1[m-2]..t1[m])+

e*(Int(p(x)*diff(func_piece1(x,m-1),x)*diff(func_piece1(x,m),x),x=t1[m-1]..t1[m])
+Int(q(x)*func_piece1(x,m-1)*func_piece1(x,m),x=t1[m-1]..t1[m]))):

for k from 2 by 1 to m-1
do

r1[k]:=evalf(

c1[k-1]*(Int(p(x)*diff(func_piece1(x,k-1),x)*diff(func_piece1(x,k),x),x=t1[k-1]..t1[k]) +Int(q(x)*func_piece1(x,k-1)*func_piece1(x,k),x=t1[k-1]..t1[k]))+

c1[k]*(Int(p(x)*diff(func_piece1(x,k),x)^2,x=t1[k-1]..t1[k+1])
+Int(q(x)*func_piece1(x,k)^2,x=t1[k-1]..t1[k+1]))+

c1[k+1]*(Int(p(x)*diff(func_piece1(x,k+1),x)*diff(func_piece1(x,k),x),x=t1[k]..t1[k+1])
+Int(q(x)*func_piece1(x,k)*func_piece1(x,k+1),x=t1[k]..t1[k+1]))-
Int(func_piece1(x,k)*f(x),x=t1[k-1]..t1[k+1]))

end do:

for i from 1 to m
do
c1[i]:=evalf( solve(r1[i], c1[i]), 15):
end do:

y_20(x):=c1[0]*func_piece1(x,0)+sum(c1[j]*func_piece1(x,j),j=1..m):

`>`	#С помощью команды dsolve получим численное решение исследуемого уравнения для сравнения с полученным по методу конечных элементов.

`>`	ds := dsolve({-diff(diff(y(x),x)p(x),x)+q(x)y(x)=f(x), y(a)=d, y(b)=e}, numeric, output=listprocedure):

`>`	#В данном цикле вычислим значения аналитической функции в узлах x[i].

`>`	for i from 0 to m do x[i] := a+(b-a)/m*i: p := eval(y(x),ds): evalf(p,5) end do:

`>`	#Полученные вычисления представим в виде таблицы.

`>`	delta_16:=0:delta_21:=0:

         for i from 1 to m
                            do
                              x[i]:=a+(b-a)/m*i;
       evalf(print([evalf(x[i],3),evalf(p(x[i],5)),simplify(subs(x=x[i],y_10(x))),
                                     simplify(subs(x=x[i],y_20(x)))]),5):
       delta_16:= delta_16+( p(x[i])-simplify(subs( x=x[i], y_10(x) )) )^2:
       delta_21:= delta_21+( p(x[i])-simplify(subs( x=x[i], y_20(x)) ) )^2:

end do:

`>`	#Командой odeplot зададим аналитическое решение.

`>`	odeplot(ds,color=blue, style=POINT, thickness=5, frames=100,view=[3..7,-0.01..1.62]);

`>`	#Командой print выведем результаты на экран.

`>`	print(sqrt(delta_16),sqrt(delta_21));

`>`	#Используем команду animate для анимирования результатов.

`>`	animate(plot,[[y_10(x), y_20(x)], x=3..t,color=[green,red],style=point], t=3..7,frames=150, view=[3..7,-0.01..1.62]);

[Maple Plot]

`>`	#Анимирование результатов вычислений для функций 10 и 20 порядка

`>`	p:=animate(plot,[[y_10(x), y_20(x)], x=3..t,color=[green,red],style=point,thickness=[3,3]],t=3..7,frames=150, view=[3..7, -0.1..1.62],scaling=constrained, background=odeplot(ds,color=black)): delta := 0.05:

`>`	#Функцией textplot выведем на экран текстовые подписи графиков.

`>`	t1 := textplot([4.3,0.8,`numeric`],color=black): t2 := textplot([4.8,0.8,`y_10(x)`],color=green): t3 := textplot([5.3,0.8,`y_20(x)`],color=red): #Функцией display выведем на экран результаты вычислений p и текстовую часть display({p,t1,t2,t3}).