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

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

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

ЧГУ, матфак, 4 курс Осипов.П.Г

                                                                          www.domath.ru

head_programm

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

>    restart:
with(plots):
 
>    #Входные данные [a,b]-отрезок, с - значение функции при x=a

>    a:=0:b:=1:
         c:=1:

#Зададим известные функции
           f:=x -> exp(2*x):
           p:=x -> -3:
#Аналитическое решение уравнения dy/dx=p(x)y+f(x)
#Такое решение может быть полученно с помощью функции dsolve, однако, я не поленился и решил
#такое уравнение в ручную
.
                     y_x:=x -> 1/5*exp(2*x)+4/5*exp(-3*x):  
#real function

body_program

>    #Процедура от двух параметров, x - переменная, n - число узлов (целое число)

>    y:=proc(x,n::integer)
                  local i,t::integer,func_piece,k,c,r; c[0]:=1:
#Разобьем отрезок на n частей
for i from 0 to n
                    do
                        t[i]:=a+i*(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;
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    
#Зададим систему алгебраических уравнений
for i from 1 to n
                 do


  for k from 1 to n-1
                     do
                       r[k]:=
c[k-1]*int( (diff(func_piece(x,k-1),x)-p(x)*func_piece(x,k-1))*func_piece(x,k),x=t[k-1]..t[k])
+
c[k]*int( (diff(func_piece(x,k),x)-p(x)*func_piece(x,k) )*func_piece(x,k), x=t[k-1]..t[k+1])
+
c[k+1]*int( (diff(func_piece(x,k+1),x)-p(x)*func_piece(x,k+1) )*func_piece(x,k), x=t[k]..t[k+1])-int(f(x)*func_piece(x,k),x=t[k-1]..t[k+1]);

r[n]:=
c[n-1]*int( (diff(func_piece(x,n-1),x)-p(x)*func_piece(x,n-1) )*func_piece(x,n),x=t[n-1]..t[n])
+
c[n]*int( (diff(func_piece(x,n),x)-p(x)*func_piece(x,n) )*func_piece(x,n), x=t[n-1]..t[n])
-int(f(x)*func_piece(x,n),x=t[n-1]..t[n]);

                      end do;

#Решение системы с помощью функции Solve
#Полученную систему можно разрешить отдельно любым из известных методов.
#Решения СЛАУ можно добиться, например, методом прогонки.

                     c[i]:=evalf( solve(r[i], c[i]), 8 ):
               end do:
#Подстановка полученных вычислений в уравнение
   y(x,n):=func_piece(x,0)+sum(c[j]*func_piece(x,j),j=1..n);

#Конец процедуры
end proc:

dates_of_function

>    #Ниже заложен механизм сравнения полученных вычислений

>    delta_6:=0:delta_11:=0:#Нулевые значения

>    for i from 0 by 0.2 to 1
                        do
         print (i ,y_x(i), simplify(subs(x=i,y(x,5))),simplify(subs(x=i,y(x,20))) );
         delta_6:= delta_6+( y_x(i)-simplify(subs( x=i, y(x,5) )) )^2:
         delta_11:= delta_11+( y_x(i)-simplify(subs( x=i, y(x,20)) ) )^2:
                                                                              end do:

0, 1, 1., 1.

.2, .7374142485, .7482388600, .7375797862

.4, .6860635551, .6717813602, .6863139399

.6, .7962624952, .8027628997, .7965315104

.8, 1.063180847, 1.034219627, 1.063394669

1.0, 1.517640875, 1.518837600, 1.517697800

graphics

>    #Вывод на экран погрешностей

>    print(sqrt(delta_6),sqrt(delta_11));   

>    #Функция plot выводит на экран графики полученных решений

>    plot( [y_x(x), y(x,5), y(x,20) ], x=0..1,0.5..1.5, color=[red,blue,yellow],                                                           style=[point,line,line],
thickness=[5,5,5],legend=["y(x)","y(x,5)","y(x,20)"]);
 

.3469283757e-1, .4598094815e-3

[Maple Plot]