Тука сложи заглавие
страницата се нуждае от дописване/преглеждане
Квадратурни форми
Миналата година се запознахме с три такива - на трапеци, правоъгълници и Симпсън - за пресмятане на определения
интеграл I(f) = \Int_a^b f(x)dx
Q^{пр}(f) = (b-a)f(\frac {a + b} 2)
Q^{тр}(f) = \frac {b - a} 2 (f(a) + f(b))
Q^с(f) = \frac {b - a} 6 (f(a) + 4f(\frac {a + b} 2) + f(b))
Ако функцията притежава производни до достатъчен ред, означаваме остатъка с
R(Q;f) = I(f) - Q(f)
R(Q^{пр}; f) = \frac {f''(\xi)} {24} (b - a)^3
R(Q^{тр}; f) = -\frac {f''(\xi)} {12} (b - a)^3
R(Q^с; f) = -\frac {f^(4)(\xi)} {2880} (b - a)^5
\xi \in [a,b] - различна за различните формули и функции
Вместо тези формули, се използват съставни квадратурни формули: делим интервала на подинтервали, прилагаме
формулите и сумираме. Съответните формули са:
|R(Q_m^{пр}; f)| \le \frac M {24m^2} (b - a)^3
|R(Q_m^{тр}; f)| \le -\frac M {12m^2} (b - a)^3
|R(Q_m^с; f)| \le -\frac M {2880m^2} (b - a)^5
Типична задача за тези формули е за кое m формулата дава най-точно приближение за интеграл от дадена
функция
Задача: Да се посочат такива m, за които m-тата съставна кваратурна формула (една от трите) приближава
a) \Int_0^1 \frac 1 {1 + x} dx (= ln 2)
b) \Int_0^1 \frac 1 {1 + x^2} dx (= \frac {\pi} 4)
с точност \epsilon = 10^{-6}
Решение
a) f(x) = \frac 1 {1 + x}, f'(x) = -\frac 1 {(1+x)^2},
f''(x) = \frac 2 {(1 + x)^3}, f^{(4)} = \frac {4!} {(1 + x)^5}
\max_{x \in [0,1]} |f^(k)(x)| = k! k = 1, 2, 3, …
|R_m(Q^{пр};f)| \le \frac {2!} {24m^2} (1 - 0)^3 = \frac 1 {12m^2} \le \frac 1 {1000000}
\Leftrightarrow m^2 \gt \frac {1000000} {12} \Leftrightarrow m \ge \frac {1000} {sqrt(12)}
|R_m(Q^{тр};f)| \le \frac {2!} {12m^2} (1 - 0)^3 = \frac 1 {6m^2} \le \frac 1 {1000000}
\Leftrightarrow m \ge \frac {1000} {sqrt(6)}
|R_m(Q^с;f)| \le \frac {4!} {2880m^4} (1 - 0)^5 = \frac 1 {120m^4} \le \frac 1 {1000000}
\Leftrightarrow m^4 \ge \frac {100000} {12} m ~ 9.55 \rightarrow 10 е достатъчно
Имаме качествена разлика между броя на интервалите при правоъгълниците и трапеците, от една страна, и Симпсън,
от друга. Причината е в скоростта на клонене към нула в двата случая, но и от растежа на втората производна.
Получените стойности са само достатъчни - можем и с доста по-малки стойности на m да получим стойността на интеграла
с исканата от нас точност. В конкретния случай можем да използваме, че когато функцията има постоянна втора
производна, грешката при правоъгълниците е положителна, а при трапеците - отрицателна. Точната стойност на интеграла
се намира между резултатите от двете квадратурни формули.
Q_m^{пр}(f) < \Int_0^1 \frac 1 {1 + x} dx < Q_m^{тр} (f)
В този случай, R(Q_m^{тр, пр}) < Q_m^{тр}(f) - Q_m^{пр}(f)
За интервала [0, 1]: x_0 = 0…..x_1…..x_2….. (…) …..1 = x_m
x_k = \frac k m, k = 0, 1, …, m \cfrac {x_k + x_{k+1}} 2 = \frac {2k + 1} {2m}
Q_m^{пр}(f) = \frac 1 m \sum_{k = 0}^{m - 1} f(\frac {2k + 1} {2m})
Q_m^{тр}(f) = \frac 1 m (\frac 1 2 f(0) + \frac 1 2 f(1) + \sum_{k = 0}^{m - 1} f(\frac k m))
Q1[f_,m_] := Sum[f[(2k+1)/(2m)], {k, 0, m-1}] / m;
Q2[f_, m_] := (f[0] + f[1])/(2m) + Sum[f[k/m], {k, 1, m-1}] / m ;
f[x_] := 1/(1+x);
Do[Print[{N[Q1[f, m], 10], N[Q2[f, m], 10], N[Q2[f, m] - Q1[f, m], 10]}], {m, 5, 25}]
Do[Print[{N[Q1[f, m], 10], N[Q2[f, m], 10], N[Q2[f, m] - Q1[f, m], 10]}], {m, 50, 75}]
Do[Print[{N[Q1[f, m], 10], N[Q2[f, m], 10], N[Q2[f, m] - Q1[f, m], 10]}], {m, 200, 210}]
Do[Print[{N[Log[2] - Q1[f, m]], N[Q2[f, m] - Log[2]]}], {m, 180, 200}]
(* Получихме желаната точност (порядък 10^(-7) с m, значително по-малко от полученото по-горе *)
Сега, ще реализираме квадратурната формула на Симпсън. При него, във всяко интервалче участват крайните точки и
средата му. Вътрешните точки се срещат 2 пъти - първо като десен, после като ляв край на интервал.
x_k = \frac k {2m}, k = 0, 1, 2, …
Q_m^c (f) = \frac {b - a} {6m} (f(x_0) + f(x_{2m}) + 2\sum_{k=1}^{m-1} f(x_2k) + 4\sum_{k=0}^{m-1} f(x_{2k-1})
???
(* Добави към горния код *)
Q3[f_, m_] := (f[0] + f[1])/(6m) + Sum[f[k/m], {k, 1, m-1}] / (3m) + 2Sum[f[(2k + 1)/(2m)], {k, 0, m-1}] / (3m);
Do[Print[N[Q3[f, m] - Log[2], 15]], {m, 3, 30}]
За b) : Как да пресмятаме (\frac 1 {1 + x^2})^{(k)} ?
\frac 1 {x^2 + 1} = \frac 1 {(x - i)(x + i)} = \frac 1 {2i} (\frac 1 {x - i} - \frac 1 {x + i}) =
- \frac i 2 (\frac 1 {x - i} - \frac 1 {x + i})
(\frac 1 {1 + x^2})^{(k)} = - \frac i 2 (\cfrac {(-1)^k k!} {(x - i)^{k+1}} - \cfrac {(-1)^k k!} {(x + i) ^{k+1}}) =
\cfrac {(-1)^{k+1} k! i} {2(x+1)^{k+1}} ((x+i)^{k+1} - (x-i_^{k+1})
(\frac 1 {x^2 + 1})'' = \frac i {(x^2+1)^3} ((x+1)^3 - (x-1)^3)
(x + i)^3 = x^3 + 3x^2i - 3x -i
(x - i)^3 = x^3 - 3x^2i - 3x +i
(x + i)^3 - (x - i)^3 = 6x^2i - 2i
\cfrac {i.i(6x^2 - 2} {(x^2 + 1)^3} = \cfrac {2(1 - 3x^2)} {{(x^2 + 1)^3}
\max_{x \in [0,1]} |f''(x)| = |f''(0)| = 2
Резултатите от a) а правоъгълниците и трапеците са същите
Гаусови квадратурни формули
Съдържанието на лекцията е изпратено по пощата преди празниците.
$Q(f) = \sum_{x=1}^n a_k f(x_k)$ - а приближено пресмятане на $I(f) = \int_a^b \mu(x) f(x) dx$
Q е Гаукова (има АСТ = 2n-1 (алгебрическа степен на точност))
${x_k}_1^n$ - нули на n-тия ортогонален полином в [a,b] при тези $\mu(x)$.
Задача. Да се намери Гаусовата формула с 2 възела за интервала [0,1] при тегло $\mu(x) \equiv 1$ .
$Q(f) = a_1 f(x_1) + a_2 f(x_2)$ , $x_1$, $x_2$ - нули на $P(x) = x^2 + ax + b$
$\int_0^1 (x^2 + ax +b) . 1 dx = 0$ //ортогонален на полиномите от степен 0
$\int_0^1 (x^2 + ax +b) . x dx = 0$ //ортогонален на полиномите от степен 1
->
$\frac 1 3 + \frac 1 2 a + b = 0$
$\frac 1 4 + \frac 1 3 a + \frac 1 2 b = 0$
->
$\frac 1 3 - \frac 1 2 + \frac 1 2 a - \frac 2 3 a = 0$
$- \frac 1 6 - \frac 1 6 a = 0 , a = -1$
$b + \frac 1 3 - \frac 1 2 = 0 , b = \frac 1 6$
$P(x) = x^2 - x + \frac 1 6 = 0, x_{1,2} = \frac {1 \plusminus sqrt(\frac 1 3)} 2$
$Q(f) = a_1 f(\frac {3 - \sqrt(3)} 6) + a_2 f(\frac {3 + \sqrt(3)} 6)$
$a_1 a_2$ намираме от условието $Q(f) = I(f)$ за f(x) = 1, x .
$f(x) = 1$ , $Q(f) = a_1.1 + a_2.1 = \int_0^1 1.dx = 1 , a_1 + a_2 = 1$
$f(x) = x$ , $Q(f) = a_1 \frac {3 - sqrt(3)} 6 + a_2 \frac {3 + sqrt(3)} 6 = \int_0^1 xdx = \frac 1 2$ ,
$a_1 {\frac (3 - sqrt(3))} 6 + a_2 {\frac (3 + sqrt(3))} 6 = \frac 1 2$
Получава се a$_1 = a_2 =\frac 1 2$
$Q(f) = \frac 1 2 (f(\frac {3 - \sqrt(3)} 6) + f(\frac {3 + \sqrt(3)} 6)) , АСТ(Q) = 3$
Квадратурни формули от Гаусов тип
За I(f) = \int_a^b \mu(x) f(x) dx
Лява квадратурна формула на Радо
Q(f) = B.f(a) + \sum_{k=1}^n a_k f(x_k) АСТ(Q) = 2
{x_k}_1^n - нули на n-тия ортогонален полином в [a, b] при тегло \mu(x) (x - a)
Дясна квадратурна формула на Радо
Q(f) = B.f(b) + \sum_{k=1}^n a_k f(x_k) АСТ(Q) = 2
{x_k}_1^n - нули на n-тия ортогонален полином в [a, b] при тегло \mu(x) (b - x)
Квадратурна формула на Лобато
Фиксирани два края на интервала на интегриране.
$Q(f) = B_1 f(a) + B_2 f(b) + \sum_{k=1}^n a_k f(x_k) АСТ(Q) = 2n + 1$
${x_k}_1^n$ - Нули на n-тия ортогонален полином в [a, b] при тегло $\mu(x)(x - a)(b - x)$