#` Задаем функцию веса p(x) и ее пределы интегрирования`
p := ln(x);  a := 0;  b := 1/2;
#` Вычисляем моменты функции p(x)`
M0 := int(p, x = a .. b);  M1 := int(p * x, x = a .. b);
M2 := int(p * x^2, x = a .. b);  M3 := int(p * x^3, x = a .. b);
#` Записываем систему уравнений на коэффициенты s и r`
eq1 := M1 * r + M2 * s = -M0;  eq2 := M2 * r + M3 * s = -M1;
#` Решаем систему уравнений для s и r`
sol := solve({eq1, eq2}, {s, r}) assuming r > 0 and s > 0;
#` Находим корни многочлена omega(x) = (x - x1) (x - x2)`
x1, x2 := solve(x^2 + sol[r] * x + sol[s], [x]);
#` Определяем коэффициенты квадратурной формулы`
A1, A2 := solve({x1 + x2 = b - a, x1 * x2 = -sol[r]}, {A1, A2});
#` Выводим результаты`  print("M0 = ", evalf(M0));  print("M1 = ",
evalf(M1));
 print("M2 = ", evalf(M2));
print("M3 = ", evalf(M3));
print("s = ", evalf(sol[s]));
  print("r = ", evalf(sol[r]));
print("x1 = ", evalf(x1));
print("x2 = ", evalf(x2));
print("A1 = ", evalf(A1));
print("A2 = ", evalf(A2));
 #` Проверяем точность квадратурной формулы на многочленах степени не выше 3`
error := int(p * x^3, x = a .. b) - evalf(A1 * (x1 - a)^3 * p(x1) + A2 * (x2 - a)^3 * p(x2));
print("Error for polynomials of degree 3 or less: ",
simplify(error));