Квадратурные оценочные точки и веса Гаусса-Лагера

Открыть Live Script

Этот пример показывает, как решить полиномиальные уравнения и системы уравнений и работать с результатами с помощью Symbolic Math Toolbox™.

Гауссовы квадратурные правила аппроксимируют интеграл по суммам $\int_{a}^{b} f (t) w (t) d t \approx \sum_{i = 1}^{n} f (x_{i}) α_{i}$ . Здесь, $x_{i}$ и $α_{i}$ являются параметрами метода, в зависимости от $n$ но не на $f$ . Они следуют из выбора функции веса $w (t)$ , следующим образом. Связанная с весовой функцией является семейством ортогональных полиномов. Корнями полиномов являются точки оценки $x_{i}$ . Наконец, веса $α_{i}$ определяются условием, что метод будет правильным для полиномов малой степени. Рассмотрим функцию веса $w (t) = \exp (- t)$ на интервале $[0, \infty]$ . Этот случай известен как квадратура Гаусса-Лагера.

syms t
n = 4;
w(t) = exp(-t);

Предположим, что вы знаете первое $n$ представители семейства ортогональных полиномов. В случае рассматриваемого здесь квадратурного правила они оказываются полиномами Лагера.

F = laguerreL(0:n-1, t)

F = 
 $(\begin{array}{c} 1 & 1 - t & \frac{t^{2}}{2} - 2 t + 1 & - \frac{t^{3}}{6} + \frac{3 t^{2}}{2} - 3 t + 1 \end{array})$

Позвольте L быть $n + 1$ st полинома, коэффициенты которого еще предстоит определить.

X = sym('X', [1, n+1])

X =  $(\begin{array}{c} X_{1} & X_{2} & X_{3} & X_{4} & X_{5} \end{array})$

L = poly2sym(X, t)

L =  $X_{1} t^{4} + X_{2} t^{3} + X_{3} t^{2} + X_{4} t + X_{5}$

Представление отношений ортогональности между полиномами Лагера F и L в системе уравнений sys.

sys = [int(F.*L.*w(t), t, 0, inf) == 0]

sys =  $(\begin{array}{c} 24 X_{1} + 6 X_{2} + 2 X_{3} + X_{4} + X_{5} = 0 & - 96 X_{1} - 18 X_{2} - 4 X_{3} - X_{4} = 0 & 144 X_{1} + 18 X_{2} + 2 X_{3} = 0 & - 96 X_{1} - 6 X_{2} = 0 \end{array})$

Добавьте условие, что полином имеет норму 1.

sys = [sys, int(L^2.*w(t), 0, inf) == 1]

sys =  $(\begin{array}{c} 24 X_{1} + 6 X_{2} + 2 X_{3} + X_{4} + X_{5} = 0 & - 96 X_{1} - 18 X_{2} - 4 X_{3} - X_{4} = 0 & 144 X_{1} + 18 X_{2} + 2 X_{3} = 0 & - 96 X_{1} - 6 X_{2} = 0 & 40320 {X_{1}}^{2} + 10080 X_{1} X_{2} + 1440 X_{1} X_{3} + 240 X_{1} X_{4} + 48 X_{1} X_{5} + 720 {X_{2}}^{2} + 240 X_{2} X_{3} + 48 X_{2} X_{4} + 12 X_{2} X_{5} + 24 {X_{3}}^{2} + 12 X_{3} X_{4} + 4 X_{3} X_{5} + 2 {X_{4}}^{2} + 2 X_{4} X_{5} + {X_{5}}^{2} = 1 \end{array})$

Решить для коэффициентов L.

S = solve(sys, X)

S = struct with fields:
    X1: [2x1 sym]
    X2: [2x1 sym]
    X3: [2x1 sym]
    X4: [2x1 sym]
    X5: [2x1 sym]

solve возвращает два решения в массиве структур. Отобразите решения.

structfun(@display, S)

ans = 
 $(\begin{array}{c} - \frac{1}{24} \\ \frac{1}{24} \end{array})$

ans = 
 $(\begin{array}{c} \frac{2}{3} \\ - \frac{2}{3} \end{array})$

ans = 
 $(\begin{array}{c} - 3 \\ 3 \end{array})$

ans = 
 $(\begin{array}{c} 4 \\ - 4 \end{array})$

ans = 
 $(\begin{array}{c} - 1 \\ 1 \end{array})$

Сделать решение уникальным путем наложения дополнительного условия, чтобы первый коэффициент был положительным:

sys = [sys, X(1)>0];
S = solve(sys, X)

S = struct with fields:
    X1: [1x1 sym]
    X2: [1x1 sym]
    X3: [1x1 sym]
    X4: [1x1 sym]
    X5: [1x1 sym]

Замените решение на L.

L = subs(L, S)

L = 
 $\frac{t^{4}}{24} - \frac{2 t^{3}}{3} + 3 t^{2} - 4 t + 1$

Как ожидалось, этот полином является |n'th многочленом Лагера:

laguerreL(n, t)

ans = 
 $\frac{t^{4}}{24} - \frac{2 t^{3}}{3} + 3 t^{2} - 4 t + 1$

Точки оценки $x_{i}$ являются корнями полинома L. Решите L для точек оценки. Корни выражены в терминах root функция.

x = solve(L)

x = 
 $\begin{array}{l} (\begin{array}{c} root (σ_{1}, z, 1) \\ root (σ_{1}, z, 2) \\ root (σ_{1}, z, 3) \\ root (σ_{1}, z, 4) \end{array}) \\ where \\ σ_{1} = z^{4} - 16 z^{3} + 72 z^{2} - 96 z + 24 \end{array}$

Форма решений может свидетельствовать о том, что ничего не достигнуто, однако на них имеются различные операции. Вычислите приближения с плавающей точкой, используя vpa:

vpa(x)

ans = 
 $(\begin{array}{c} 0.32254768961939231180036145910437 \\ 1.7457611011583465756868167125179 \\ 4.5366202969211279832792853849571 \\ 9.3950709123011331292335364434205 \end{array})$

Могут возникнуть паразитные мнимые части. Докажи символически, что корни являются вещественными числами:

isAlways(in(x, 'real'))

ans = 4x1 logical array

   1
   1
   1
   1

Для полиномов степени, меньшей или равной 4, можно использовать MaxDegree получить решения с точки зрения вложенных радикалов вместо этого с точки зрения root. Однако последующие операции по результатам этой формы будут происходить медленно.

xradical = solve(L, 'MaxDegree', 4)

xradical = 
 $\begin{array}{l} (\begin{array}{c} 4 - σ_{1} - σ_{3} \\ 4 + σ_{1} - σ_{3} \\ σ_{3} - σ_{2} + 4 \\ σ_{3} + σ_{2} + 4 \end{array}) \\ where \\ σ_{1} = \frac{\sqrt{96 σ_{6} σ_{4} - 3 σ_{5} σ_{4} - 288 σ_{4} - 512 \sqrt{3} \sqrt{6} \sqrt{6 + \sqrt{3} \sqrt{6} i}}}{2 {(768 + 128 \sqrt{3} \sqrt{6} i)}^{1 / 6} {(144 σ_{6} + 9 σ_{5} + 864)}^{1 / 4}} \\ σ_{2} = \frac{\sqrt{96 σ_{6} σ_{4} - 3 σ_{5} σ_{4} - 288 σ_{4} + 512 \sqrt{3} \sqrt{6} \sqrt{6 + \sqrt{3} \sqrt{6} i}}}{2 {(768 + 128 \sqrt{3} \sqrt{6} i)}^{1 / 6} {(144 σ_{6} + 9 σ_{5} + 864)}^{1 / 4}} \\ σ_{3} = \frac{σ_{4}}{2 {(768 + 128 \sqrt{3} \sqrt{6} i)}^{1 / 6}} \\ σ_{4} = \sqrt{16 σ_{6} + σ_{5} + 96} \\ σ_{5} = {(768 + 128 \sqrt{3} \sqrt{6} i)}^{2 / 3} \\ σ_{6} = {(768 + 128 \sqrt{3} \sqrt{6} i)}^{1 / 3} \end{array}$

Веса $α_{i}$ заданы условием, что для полиномов степени меньше, чем $n$ квадратурное правило должно привести к точным результатам. Достаточно, если это сохраняет для базиса векторного пространства этих полиномов. Это условие приводит к системе четырех уравнений в четырех переменных.

y = sym('y', [n, 1]);
sys = sym(zeros(n));
for k=0:n-1 
    sys(k+1) = sum(y.*(x.^k)) == int(t^k * w(t), t, 0, inf); 
end
sys

sys = 
 $\begin{array}{l} (\begin{array}{c} y_{1} + y_{2} + y_{3} + y_{4} = 1 & 0 & 0 & 0 \\ y_{1} root (σ_{1}, z, 1) + y_{2} root (σ_{1}, z, 2) + y_{3} root (σ_{1}, z, 3) + y_{4} root (σ_{1}, z, 4) = 1 & 0 & 0 & 0 \\ y_{1} {root (σ_{1}, z, 1)}^{2} + y_{2} {root (σ_{1}, z, 2)}^{2} + y_{3} {root (σ_{1}, z, 3)}^{2} + y_{4} {root (σ_{1}, z, 4)}^{2} = 2 & 0 & 0 & 0 \\ y_{1} {root (σ_{1}, z, 1)}^{3} + y_{2} {root (σ_{1}, z, 2)}^{3} + y_{3} {root (σ_{1}, z, 3)}^{3} + y_{4} {root (σ_{1}, z, 4)}^{3} = 6 & 0 & 0 & 0 \end{array}) \\ where \\ σ_{1} = z^{4} - 16 z^{3} + 72 z^{2} - 96 z + 24 \end{array}$

Решить систему как численно, так и символически. Решение является желаемым вектором весов $α_{i}$ .

[a1, a2, a3, a4] = vpasolve(sys, y)

a1 =  $0.60315410434163360163596602381808$

a2 =  $0.35741869243779968664149201745809$

a3 =  $0.03888790851500538427243816815621$

a4 =  $0.00053929470556132745010379056762059$

[alpha1, alpha2, alpha3, alpha4] = solve(sys, y)

alpha1 = 
 $\begin{array}{l} - \frac{σ_{3} σ_{2} + σ_{3} σ_{1} + σ_{2} σ_{1} - σ_{3} σ_{2} σ_{1} - 2 σ_{3} - 2 σ_{2} - 2 σ_{1} + 6}{(σ_{4} - σ_{3}) (σ_{4} σ_{2} + σ_{4} σ_{1} - σ_{2} σ_{1} - {σ_{4}}^{2})} \\ where \\ σ_{1} = root (z^{4} - 16 z^{3} + 72 z^{2} - 96 z + 24, z, 4) \\ σ_{2} = root (z^{4} - 16 z^{3} + 72 z^{2} - 96 z + 24, z, 3) \\ σ_{3} = root (z^{4} - 16 z^{3} + 72 z^{2} - 96 z + 24, z, 2) \\ σ_{4} = root (z^{4} - 16 z^{3} + 72 z^{2} - 96 z + 24, z, 1) \end{array}$

alpha2 = 
 $\begin{array}{l} \frac{root (σ_{1}, z, 1) root (σ_{1}, z, 3) + root (σ_{1}, z, 1) root (σ_{1}, z, 4) + root (σ_{1}, z, 3) root (σ_{1}, z, 4) - root (σ_{1}, z, 1) root (σ_{1}, z, 3) root (σ_{1}, z, 4) - 2 root (σ_{1}, z, 1) - 2 root (σ_{1}, z, 3) - 2 root (σ_{1}, z, 4) + 6}{(root (σ_{1}, z, 2) - root (σ_{1}, z, 1)) (root (σ_{1}, z, 2) - root (σ_{1}, z, 3)) (root (σ_{1}, z, 2) - root (σ_{1}, z, 4))} \\ where \\ σ_{1} = z^{4} - 16 z^{3} + 72 z^{2} - 96 z + 24 \end{array}$

alpha3 = 
 $\begin{array}{l} \frac{σ_{3} σ_{2} + σ_{3} σ_{1} + σ_{2} σ_{1} - σ_{3} σ_{2} σ_{1} - 2 σ_{3} - 2 σ_{2} - 2 σ_{1} + 6}{(σ_{4} - σ_{1}) (σ_{3} σ_{2} - σ_{3} σ_{4} - σ_{2} σ_{4} + {σ_{4}}^{2})} \\ where \\ σ_{1} = root (z^{4} - 16 z^{3} + 72 z^{2} - 96 z + 24, z, 4) \\ σ_{2} = root (z^{4} - 16 z^{3} + 72 z^{2} - 96 z + 24, z, 2) \\ σ_{3} = root (z^{4} - 16 z^{3} + 72 z^{2} - 96 z + 24, z, 1) \\ σ_{4} = root (z^{4} - 16 z^{3} + 72 z^{2} - 96 z + 24, z, 3) \end{array}$

alpha4 = 
 $\begin{array}{l} - \frac{σ_{3} σ_{2} + σ_{3} σ_{1} + σ_{2} σ_{1} - σ_{3} σ_{2} σ_{1} - 2 σ_{3} - 2 σ_{2} - 2 σ_{1} + 6}{{σ_{4}}^{2} σ_{3} + {σ_{4}}^{2} σ_{2} + {σ_{4}}^{2} σ_{1} - {σ_{4}}^{3} + σ_{3} σ_{2} σ_{1} - σ_{3} σ_{2} σ_{4} - σ_{3} σ_{1} σ_{4} - σ_{2} σ_{1} σ_{4}} \\ where \\ σ_{1} = root (z^{4} - 16 z^{3} + 72 z^{2} - 96 z + 24, z, 3) \\ σ_{2} = root (z^{4} - 16 z^{3} + 72 z^{2} - 96 z + 24, z, 2) \\ σ_{3} = root (z^{4} - 16 z^{3} + 72 z^{2} - 96 z + 24, z, 1) \\ σ_{4} = root (z^{4} - 16 z^{3} + 72 z^{2} - 96 z + 24, z, 4) \end{array}$

Также можно получить решение как структуру, задав только один выходной аргумент.

S = solve(sys, y)

S = struct with fields:
    y1: [1x1 sym]
    y2: [1x1 sym]
    y3: [1x1 sym]
    y4: [1x1 sym]

structfun(@double, S)

Преобразуйте структуру S в символьный массив:

Scell = struct2cell(S);
alpha = transpose([Scell{:}])

alpha = 
 $\begin{array}{l} (\begin{array}{c} - \frac{σ_{3} + σ_{15} + σ_{2} - σ_{6} - σ_{12} - σ_{11} - σ_{10} + 6}{(root (σ_{16}, z, 1) - root (σ_{16}, z, 2)) (σ_{1} + σ_{4} - σ_{2} - {root (σ_{16}, z, 1)}^{2})} \\ \frac{σ_{1} + σ_{4} + σ_{2} - σ_{7} - σ_{13} - σ_{11} - σ_{10} + 6}{(root (σ_{16}, z, 2) - root (σ_{16}, z, 1)) (root (σ_{16}, z, 2) - root (σ_{16}, z, 3)) (root (σ_{16}, z, 2) - root (σ_{16}, z, 4))} \\ \frac{σ_{5} + σ_{4} + σ_{15} - σ_{8} - σ_{13} - σ_{12} - σ_{10} + 6}{(root (σ_{16}, z, 3) - root (σ_{16}, z, 4)) (σ_{5} - σ_{1} - σ_{3} + {root (σ_{16}, z, 3)}^{2})} \\ - \frac{σ_{5} + σ_{1} + σ_{3} - σ_{9} - σ_{13} - σ_{12} - σ_{11} + 6}{σ_{14} root (σ_{16}, z, 1) + σ_{14} root (σ_{16}, z, 2) + σ_{14} root (σ_{16}, z, 3) - {root (σ_{16}, z, 4)}^{3} + σ_{9} - σ_{8} - σ_{7} - σ_{6}} \end{array}) \\ where \\ σ_{1} = root (σ_{16}, z, 1) root (σ_{16}, z, 3) \\ σ_{2} = root (σ_{16}, z, 3) root (σ_{16}, z, 4) \\ σ_{3} = root (σ_{16}, z, 2) root (σ_{16}, z, 3) \\ σ_{4} = root (σ_{16}, z, 1) root (σ_{16}, z, 4) \\ σ_{5} = root (σ_{16}, z, 1) root (σ_{16}, z, 2) \\ σ_{6} = root (σ_{16}, z, 2) root (σ_{16}, z, 3) root (σ_{16}, z, 4) \\ σ_{7} = root (σ_{16}, z, 1) root (σ_{16}, z, 3) root (σ_{16}, z, 4) \\ σ_{8} = root (σ_{16}, z, 1) root (σ_{16}, z, 2) root (σ_{16}, z, 4) \\ σ_{9} = root (σ_{16}, z, 1) root (σ_{16}, z, 2) root (σ_{16}, z, 3) \\ σ_{10} = 2 root (σ_{16}, z, 4) \\ σ_{11} = 2 root (σ_{16}, z, 3) \\ σ_{12} = 2 root (σ_{16}, z, 2) \\ σ_{13} = 2 root (σ_{16}, z, 1) \\ σ_{14} = {root (σ_{16}, z, 4)}^{2} \\ σ_{15} = root (σ_{16}, z, 2) root (σ_{16}, z, 4) \\ σ_{16} = z^{4} - 16 z^{3} + 72 z^{2} - 96 z + 24 \end{array}$

Символическое решение выглядит сложным. Упростите его и преобразуйте в вектор с плавающей точкой:

alpha = simplify(alpha)

alpha = 
 $\begin{array}{l} (\begin{array}{c} \frac{{root (σ_{1}, z, 1)}^{2}}{72} - \frac{29 root (σ_{1}, z, 1)}{144} + \frac{2}{3} \\ \frac{{root (σ_{1}, z, 2)}^{2}}{72} - \frac{29 root (σ_{1}, z, 2)}{144} + \frac{2}{3} \\ \frac{{root (σ_{1}, z, 3)}^{2}}{72} - \frac{29 root (σ_{1}, z, 3)}{144} + \frac{2}{3} \\ \frac{{root (σ_{1}, z, 4)}^{2}}{72} - \frac{29 root (σ_{1}, z, 4)}{144} + \frac{2}{3} \end{array}) \\ where \\ σ_{1} = z^{4} - 16 z^{3} + 72 z^{2} - 96 z + 24 \end{array}$

vpa(alpha)

ans = 
 $(\begin{array}{c} 0.60315410434163360163596602381808 \\ 0.35741869243779968664149201745809 \\ 0.03888790851500538427243816815621 \\ 0.00053929470556132745010379056762059 \end{array})$

Увеличьте читаемость путем замены вхождений корней x в alpha по сокращениям:

subs(alpha, x, sym('R', [4, 1]))

ans = 
 $(\begin{array}{c} \frac{{R_{1}}^{2}}{72} - \frac{29 R_{1}}{144} + \frac{2}{3} \\ \frac{{R_{2}}^{2}}{72} - \frac{29 R_{2}}{144} + \frac{2}{3} \\ \frac{{R_{3}}^{2}}{72} - \frac{29 R_{3}}{144} + \frac{2}{3} \\ \frac{{R_{4}}^{2}}{72} - \frac{29 R_{4}}{144} + \frac{2}{3} \end{array})$

Суммируйте веса, чтобы показать, что их сумма равна 1:

simplify(sum(alpha))

ans =  $1$

Другой метод получения весов квадратурного правила - вычислить их с помощью формулы $α_{i} = \int_{a}^{b} w (t) \prod_{j \neq i} \frac{t - x_{j}}{x_{i} - x_{j}} d t$ . Сделайте это для $i = 1$ . Это приводит к тому же результату, что и другой метод:

int(w(t) * prod((t - x(2:4)) ./ (x(1) - x(2:4))), t, 0, inf)

ans = 
 $\frac{{root (z^{4} - 16 z^{3} + 72 z^{2} - 96 z + 24, z, 1)}^{2}}{72} - \frac{29 root (z^{4} - 16 z^{3} + 72 z^{2} - 96 z + 24, z, 1)}{144} + \frac{2}{3}$

Квадратурное правило дает точные результаты даже для всех полиномов со степенью, меньшей или равной $2 n - 1$ , но не для $t^{2 n}$ .

simplify(sum(alpha.*(x.^(2*n-1))) -int(t^(2*n-1)*w(t), t, 0, inf))

ans =  $0$

simplify(sum(alpha.*(x.^(2*n))) -int(t^(2*n)*w(t), t, 0, inf))

ans =  $- 576$

Примените квадратурное правило к косинусу и сравните с точным результатом:

vpa(sum(alpha.*(cos(x))))

ans =  $0.50249370546067059229918484198931$

int(cos(t)*w(t), t, 0, inf)

ans = 
 $\frac{1}{2}$

Для степеней косинуса ошибка колеблется между нечетной и четной степенями:

errors = zeros(1, 20);
for k=1:20 
    errors(k) = double(sum(alpha.*(cos(x).^k)) -int(cos(t)^k*w(t), t, 0, inf));
end
plot(real(errors))

Figure contains an axes. The axes contains an object of type line.

Документация

Квадратурные оценочные точки и веса Гаусса-Лагера

Документация Symbolic Math Toolbox

Поддержка