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:

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

>	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);