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

Открыть сценарий в реальном времени

В этом примере показано, как решать полиномиальные уравнения и системы уравнений и работать с результатами с помощью символьных математических Toolbox™.

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

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-м многочленом, коэффициенты которого ещё предстоит определить.

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$

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

laguerreL(n, t)

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

Точки оценки $_{xi}$ являются корнями многочлена 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=∫abw (t \underset{}{)} \frac{_{}}{_{}_{}}$ ∏j≠it-xjxi-xjdt. Сделайте это $для 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}$

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

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.

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

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

Документация по инструментам символьной математики

Поддержка