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

Вычислительная математика

Метод Гаусса. Вычисление интеграла

матфак, Иванова М.

                                             www.domath.ru

>    #Объявим исходную функцию

>    f:=x->x*cos(x)/(1+x*x);

#пусть известно,что весовая функция p(x) равна 1  

      p:=x->1;

  #объявим константу, равную 4

      n:=4;

>    #Обьявим процедуру, вычисляющую интеграл от выражения p(x)*a*b,
#где a и b - это входные параметры процедуры.

#здесь 0 и 1 пределы интегрирования

>    sp:=proc(a,b)

                                                 int(p(x)*a*b,x=0..1):

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

f := proc (x) options operator, arrow; x*cos(x)/(1+x*x) end proc

p := proc (x) options operator, arrow; 1 end proc

n := 4

>    #Ниже в двойном цикле заложен механизм
#задания функции

>    for i from 1 to n do                                      

             for k from 1 to n do
      
 #если окажется, что (k>i), то все параметры a[-i,k]:=0
        if (k>i) then a[-i,k]:=0 end if;
       end do;
                                                                                                                                                                                                                                                                                                                               

                     f[i](x):=x^(n-i)+a[n-1-i,i]*x^(n-1-i)+a[n-2-i,i]*x^(n-2-i)+a[n-3-i,i]*x^(n-3-i);                                                                                                                                                                                                    end do;

f[1](x) := x^3+a[2,1]*x^2+a[1,1]*x+a[0,1]

f[2](x) := x^2+a[1,2]*x+a[0,2]

f[3](x) := x+a[0,3]

f[4](x) := 1



#здесь в двойном цикле вычисляются интегралы от выражения:
#это есть всевозможные произведения функций fi и fj,sp-заданная процедура ранее

>    for i from 1 to n do   
     for j from 1 to n do   
                        g[i,j]:=sp(f[i](x),f[j](x));
      end do;
end do;

>    #механизм гаусса

>    for i from 1 to n-1 do
                  if (i=1) then
#для первого
                           s1[i] := {g[n,n-i]};
                           s2[i] := {a[i-1,n-i]};
                   else
                       if (i=2) then
                                  
  #для второго
                               s1[i] := {g[n,n-i],g[n-1,n-i]};
                               s2[i] := {a[i-1,n-i],a[i-2,n-i]};           
                         else
                             if (i=3) then

                                 s1[i] := {g[n,n-i],g[n-1,n-i],g[n-2,n-i]};
                                 s2[i] := {a[i-1,n-i],a[i-2,n-i],a[i-3,n-i]};                                 
                              end if;
                          end if;
                    end if;

            sys[i]:=solve(s1[i],s2[i]);
 
     f[n-i](x):=subs(sys[i],f[n-i](x));

end do;

sys[1] := {a[0,3] = -1/2}

f[3](x) := x-1/2

sys[2] := {a[1,2] = -1, a[0,2] = 1/6}

f[2](x) := x^2-x+1/6

sys[3] := {a[1,1] = 3/5, a[2,1] = -3/2, a[0,1] = -1/20}

f[1](x) := x^3-3/2*x^2+3/5*x-1/20

>    decide:=[solve(f[1](x))];#выведем искомое решение

decide := [1/2, 1/2+1/10*15^(1/2), 1/2-1/10*15^(1/2)]

>    for i from 1 to n-1 do
                     x[i]:=op(i,decide);
od;

x[1] := 1/2

x[2] := 1/2+1/10*15^(1/2)

x[3] := 1/2-1/10*15^(1/2)

>    for i from 1 to n-1 do
           
 #произведение выражений
            pr[i]:= mul(((x-x[j])/(x[i]-x[j])),j=1..i-1)*mul(((x-x[j])/(x[i]-x[j])), j=i+1..3);
            c[i]:=sp(p(x),pr[i]);
#все коэффициенты ci
od;

>    Integral:=sum('c[k]*f(x[k])','k'=1..3);
#интеграл методом гаусса

>    evalf(%);#ответ

.2738182064