PA081: Programování numerických výpočtů
3. Nelineární rovnice o jedné neznámé praktické příklady
Aleš Křenek
jaro 2011
Montáž ložiska
K sestavení kluzného ložiska uložení ocelové mostní konstrukce je třeba na místě zasunout čep o vnějším průměru 313 mm do náboje o stejném vnitřním průměru.
Sestavení probíhá s podchlazeným čepem, díky tepelné roztažnosti materiálu je menší a ložisko lze sestavit. Bezpečné zmenšení průměru je 0.3 7 mm.
Na jakou teplotu je třeba čep podchladit, předpokládáme-li teplotu prostředí 20°C?
PA081: Programování numerických
výpočtů
A. Křenek
Montáž ložiska
Cirkusové číslo
Obchodní případ
2/22
Montáž ložiska
► koeficient tepelné roztažnosti oceli a závisí na teplotě
► výpočet tedy není zcela triviální
► v rozmezí ±200°C je poměrně přesně vyjádřen funkcí
a(T) = aT2 + bT + C
kde a = -8.27 x ÍO"11,*? = 2.21 x l(r8,c = 1.17 x 1(T5
► koeficienty získány regresí (proložením křivky) experimentálních dat
► metodu regrese budeme probírat později
Montáž ložiska
► poloměr po zchlazení počítáme pro AT ->■ 0
rx = r0(l-ATa(T0))(l-ATa(T0-AT))...
= r0(l - ATa(To) - ATa(T0 - AT)
+ (AT)2a(To)a(T0-AT))...
► vedlo by na netriviální diferenciální rovnici
► protože a(T) «: 1, AT «: 1, lze nelineární členy zanedbat
► dostáváme zjednodušení
rx = r0(l - &(TX)) kde a(Tx) = a budeme řešit rovnici
•7b
a(T)dT
a(Tx)
r o -rx To
Montáž ložiska
► řešená funkce
f(Tx) = a(Tx)-
r0~rx
To b, 2
-rá + -rz + ct
To
r0~rx
f(7o3 - Tx) + ^(T02 - 7*) + c(T0 - Tx)
tx r0 r o -rx To
► potřebujeme
► separaci „správného" kořene
► zvolit vhodnou numerickou metodu
PA081: Programování numerických
výpočtů
A. Křenek
Montáž ložiska
Cirkusové číslo
Obchodní případ
5/22
Montáž ložiska
► řešená funkce
f(Tx) = a(Tx)-
r0~rx
To b, 2
-rá + -rz + ct
To
r0~rx
f(7o3 - Tx) + ^(T02 - 7*) + c(T0 - Tx)
tq -rx
t0
► potřebujeme
► separaci „správného" kořene
► zvolit vhodnou numerickou metodu
► koeficient u Tx je kladný (a bylo záporné)
► pro Tx -* oo také f(Tx) -* oo
► pro Tx -> -oo také f(Tx) -> -oo
PA081: Programování numerických
výpočtů
A. Křenek
Montáž ložiska
Cirkusové číslo
Obchodní případ
5/22
Montáž ložiska
► derivace řešené funkce
f'(Tx) = -aTx-bTx-c
* f(Tx) má jedno lokální minimum a jedno maximum
► v bodech, kde/'(T%) =0
► numericky stabilní výpočet?
► druhá přednáška ;-)
► b = -2.21 x ÍCT8, -JD = 6.60 x ÍCT8 je OK
► při nepřesném určení by mohla Newtonova metoda zabloudit, viz dále
► vypočtené extrémy
Ti = -265.544   /(Ti) = 0.000867609 T2 = 532.775    f(T2) = -0.00614507
Montáž ložiska
► dva „nesmyslné" kořeny T- < Ti a T+ > T2
► nemají fyzikální význam
► T+ - čep se při vysokých teplotách začne smršťovat
► T- - při nízkých teplotách se začne znovu roztahovat
► obojí důsledkem aplikace regresní křivky mimo rozsah hodnot, ze kterých byla určena
skutečný kořen v intervalu [Ti, T2]
► funkce je spojitá, kořen je právě jeden
► dokonalá separace
► umíme jednoduše počítat derivaci
► navíc je na (ľi, T2) derivace nenulová
► Newtonova metoda je ideální
Montáž ložiska
► numerická stabilita výpočtu f (T) a f (T)
► rozsah T řádově ve stovkách
► členy polynomu v řádech 10~5-10~2
► raději použijeme typ double
► potenciální problém Newtonovy metody f (T) ->■ 0
► nenastane díky známým vlastnostem funkce
Obchodní případ
8/22
Montáž ložiska
► Průběh výpočtu funkcí rtnewt
► požadovaná absolutní přesnost 10~3
► nemá smysl více, vstupní data (koeficienty a, b, c) jsou méně přesné
► a jak budeme při stavbě mostu chladit čep na takhle přesnou teplotu?
f(T)
f (T)
133.615
-66.6454
-89.1
-90.0546
-90.0564
-0.00263873
-0.000221398
-8.66255e-06
-1.68088e-08
-6.3964e-14
-1.31765e-05 -9.85981e-06 -9.07435e-06 -9.03911e-06 -9.03904e-06
► v následujícím kroku už bylo dosaženo požadované přesnosti
► výsledek je -90.056°C
Cirkusové číslo
Artisté připravují nové vystoupení.
Na začátku se jeden z nich zachytí rukama i nohama
v kruhové konstrukci o průměru 180 cm a další konstrukci
vykoulí na plošinu.
Je třeba, aby se na konci tohoto pohybu původně nejvyšší bod kruhové konstrukce dotýkal plošiny a byl ve stejné výšce jako na začátku.
Jaká je potřebná délka a sklon plošiny?
11/22
Cirkusové číslo
označíme
► R - poloměr kruhu
► a - úhel otočení kruhu (v radiánech)
► P - úhel naklonění plošiny
délka pohybu kruhu po plošině je Ra z naznačeného čtyřúhelníku odečteme
► a + P = tt (další dva úhly jsou pravé) ^ tanf = £
z toho vyplývá řešená rovnice
Obchodní prípad
tanJ
tt
tedy     f (P) = (7t-J6)tan
2
12/22
Cirkusové číslo
► separace kořene
► P je ostrý úhel, tj. P
(0,f)
řešení je právě jedno
z geometrické podstaty problému
formulací rovnic jsme nepřidali falešný kořen
dáme si práci s exaktním důkazem nebo to rovnou zkusíme
pro jistotu ověříme /(O) = -
/(tt/2) = 0.5707963
► numericky nebezpečná by byla až oblast j6 ->• tt
► vedla by na součin typu 0 ■ oo pohybujeme se v bezpečné vzdálenosti
► budeme předstírat, že neumíme spočítat derivaci
► ukážeme metodu sečen, Riddersovu, a Brentovu
výpočtu A. Křenek
Obchodní případ
13/22
Cirkusové číslo
Metoda sečen, požadovaná absolutní přesnost 10~4
n	beta		f(beta)	
1	+0.	.0000000	-1.	.0000000
2	+1.	.5707963	+0.	.5707963
3	+0.	.7853982	-0.	.0240323
4	+1.	.1780972	+0.	.3119657
B	+0.	.9817477	+0.	.1544612
6	+0.	.8835729	+0.	.0679638
7	+0.	.8344855	+0.	.0226702
8	+0.	.8099419	-0.	.0005027
9	+0.	.8222137	+0.	.0111281
15	+0.	.8105171	+0.	.0000445
16	+0.	.8104212	-0.	.0000467
► počínaje třetím krokem výpočet vždy ve středu intervalu
► pomalá konvergence
Obchodní případ
14/22
Cirkusové číslo
Riddersova metoda, požadovaná absolutní přesnost 10
n	beta		f(beta)	
1	+0.	.0000000	-1.	.0000000
2	+1.	.5707963	+0.	.5707963
3	+0.	.7853982	-0.	.0240323
4	+0.	.8103685	-0.	.0000968
B	+1.	.1905824	+0.	.3213144
6	+0.	.8104702	-0.	.0000001
7	+1.	.0005263	+0.	.1704016
(8)	+0.	.8104702		
4. první iterační výpočet
5. půlení intervalu mezi 2. a 4.
6. další iterace
7. půlení mezi 5. a 6.
8. výsledek, už je v toleranci
Cirkusové číslo
Brentova metoda, požadovaná absolutní přesnost 10~4
n	beta		f(beta)	
1	+0.	.0000000	-1.	.0000000
2	+1.	.5707963	+0.	.5707963
3	+1.	.0000000	+0.	.1699574
4	+0.	.7931377	-0.	.0165737
B	+0.	.8115179	+0.	.0009960
6	+0.	.8104759	+0.	.0000054
7	+0.	.8104259	-0.	.0000422
(8)	+0.	.8104759		
Obchodní prípad
► počítá primárním způsobem, nedošlo na bezpečné půlení intervalů
► srovnatelně rychlá konvergence s Riddersem
► při větší požadované přesnosti se liší ± o jeden krok
16/22
Cirkusové číslo
► všechny metody se v rámci požadované tolerance shodují
► rychlost konvergence dle očekávání
► zdvojnásobení požadované přesnosti zdvojnásobí počet kroků půlení intervalů
► Riddersova metoda naroste ze 7 na 9
► Brentova dokonce jen ze 7 na 8
► lepší výsledky než teoretická rychlost konvergence V2
PA081: Programování numerických
výpočtů
A. Křenek
Montáž ložiska
Cirkusové číslo
Obchodní případ
17/22
Cirkusové číslo
► délka pohybu po plošině je
Ra = R (n - j8) = 2.098
► zbývá spočítat část plošiny od země k bodu dotyku
tan= ?     tedy     i =-R —- = 0.3861
2        l y tan((rr - j8)/2)
► nehrozí významné numerické nepřesnosti
► celková délka plošiny je 2.484 m, sklon 46.44°
PA081: Programovaní numerických
výpočtů
A. Křenek
Montáž ložiska
Cirkusové číslo
Obchodní prípad
18/22
Obchodní případ
Americká garážová firma zahajuje montáž počítačů. Kolik jich musí v prvním roce prodat, aby byla zisková?
Obchodní případ
Americká garážová firma zahajuje montáž počítačů. Kolik jich musí v prvním roce prodat, aby byla zisková?
► uvažované náklady
► fixní jednorázové (pořízení vybavení): $20000
► fixní roční (nájem, web, ...): S15000
► variabilní (komponenty, hodinová sazba práce, ...): S62
► semivariabilní (náročnější logistika atd.) S30nL5
► předpokládané příjmy
► prodej počítačů: $1500n
► slevy (množstevní, VIP zákazníci, ...): -SlOn1-5
Obchodní případ
► řešená rovnice
f (n) = T C (n) - T S (n) = 35000 - 875n + 40ti1
► první dva členy - lineární funkce
► prochází bodem (0,35000)
► osu x protíná v 35000/875 = 40
► člen 40n15 způsobí „prohnutí" nahoru
► posune kořen z bodu 40 dál
► přidá druhý kořen pro vyšší n
► zisku dosahujeme v oblasti mezi těmito kořeny
Obchodní případ
► separace 1. kořene
► podle předchozí úvahy /(40) > 0, skutečná hodnota 10119
► zkusíme /(80) = -6378, vyhovuje
► separace 2. kořene
► předchozí /(80) = -6378
► hrubý odhad zanedbáním absolutního člene
-875w + 40wL5 = 0     tedy     n=í — )  = 478.51
► zkusíme /(500) = 44713, vyhovuje
► numerická stabilita v dané separaci kořenů
► součty/rozdíly čísel s minimálním řádovým rozdílem
► reálně nás zajímá přesnost na jednotky
► žádný potenciální problém
Obchodní případ
vypočtené kořeny 62.7 a 384.0 ► výsledky v požadované přesnosti shodné všemi metodami
chování jednotlivých metod (počty iterací)
přesnost	10	-i	10	-4
kořen	1	2	1	2
půlení	11	15	21	25
Ridders	5	9	7	11
Brent	5	9	7	10
► umělá přesnost 10 4 - zdvojnásobení počtu platných číslic
► dokládá očekávanou rychlost konvergence
lineární pro půlení intervalů \ľl pro ostatní metody
Obchodní případ
22/22