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:

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