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

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

Метод Бубнова-Галеркина

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

                                              www.domath.ru

 

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

>    restart:
with(plots):

>    #Входные данные

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

#Входные параметры функции

          f:=x -> exp(2*x):
          p:=x -> -3:
#Аналитическая функция. Такая функция может быть получена при решении соответ.
#Уравнения с помощью функции dsolve

                     y_x:=x -> 1/5*exp(2*x)+4/5*exp(-3*x):
 #real function

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

>    sp:=proc(f,a,b)
          
       int(f,x=a..b);
end proc:

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

>    y:=proc(x,n::integer)
#Локальные переменные
           local i,j,k::integer,

                  r,alfa::float; alfa[0]:=1:
#Механизм Бубнова-Галеркина    
    for j from 1 to n do

#СЛАУ

for k from 1 to n do  

      r[k]:=sp( sum(alfa[i]*( i*x^(i-1)-p(x)*x^(i) )*x^(k-1),i=1..n ),a,b)
                                       -alfa[0]*sp(p(x)*x^(k-1),a,b)-sp(f(x)*x^(k-1),a,b);
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                end do:

#Решение Слау         
           alfa[j]:= evalf( solve (r[j]=0, alfa[j] ), 7 );

         end do;
#Конец процедуры - вывод решения по методу
   alfa[0]+sum(alfa[i]*x^(i),i=1..n);
end proc:

>    #Таблица решений и погрешности

>    delta_3:=0:delta_4:=0:
for i from 0 by 0.2 to 1
                        do
          print (i ,y_x(i), subs(x=i,y(x,3)), subs(x=i,y(x,4)) );
          delta_3:=delta_3+( y_x(i)-subs( x=i, y(x,3) ) )^2:
          delta_4:=delta_4+( y_x(i)-subs( x=i, y(x,4) ) )^2:
                                                                              end do:

0, 1, 1., 1.

.2, .7374142485, .7410643417, .7369924266

.4, .6860635551, .6749046955, .6857195147

.6, .7962624952, .7897881337, .7970060045

.8, 1.063180847, 1.073981727, 1.063098524

1.0, 1.517640875, 1.515752550, 1.517665589

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

>    print(sqrt(delta_3),sqrt(delta_4));

>    #Вывод графиков решений

>    plot( [y_x(x), y(x,3), y(x,4) ], x=0..1,0.2..1.2, color=[yellow,blue,red],                                                           legend=["y(x)","y(x,3)","y(x,4)"]);y(x):=y_x(x);y(x,3):=y(x,3);y(x,4):=y(x,4);

 

.1732008848e-1, .9254682325e-3

[Maple Plot]

y(x) := 1/5*exp(2*x)+4/5*exp(-3*x)

y(x,3) := 1-1.796173202*x+2.556361752*x^2-.2444360*x^3

y(x,4) := 1-1.978110647*x+3.735999817*x^2-2.318918581*x^3+1.078695*x^4