Дисциплина:
Раздел:
Выполнил:

Методы конечных элементов

Одномерные конечные элементы. Уравнение Штурма-Лиувилля

ЧГУ, матфак, 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]);

.6367042555e-1, .1176306188e-2

[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}).