PA081: Programování numerických výpoˇ ct˚ u A. Kˇ renek Montáž ložiska
PA081: Programování numerických výpoˇ ct˚ u
Cirkusové ˇ císlo Obchodní pˇ rípad
4. Nelineární rovnice o jedné neznámé praktické pˇ ríklady
Aleš Kˇ renek
jaro 2012
1/22
Montáž ložiska
PA081: Programování numerických výpoˇ ct˚ u A. Kˇ renek Montáž ložiska Cirkusové ˇ císlo
K sestavení kluzného ložiska uložení ocelové mostní konstrukce je tˇ reba na místˇ e zasunout ˇ cep o vnˇ ejším pr˚ umˇ eru 313 mm do náboje o stejném vnitˇ rním pr˚ umˇ eru.
Obchodní pˇ rípad
Sestavení probíhá s podchlazeným ˇ cepem, díky tepelné roztažnosti materiálu je menší a ložisko lze sestavit. Bezpeˇ cné zmenšení pr˚ umˇ eru je 0.37 mm. Na jakou teplotu je tˇ reba ˇ cep podchladit, pˇ redpokládáme-li teplotu prostˇ redí 20◦ C?
2/22
PA081: Programování numerických výpoˇ ct˚ u
Montáž ložiska
A. Kˇ renek Montáž ložiska Cirkusové ˇ císlo
ñ
koeficient tepelné roztažnosti oceli α závisí na teplotˇ e ñ
ñ
výpoˇ cet tedy není zcela triviální
Obchodní pˇ rípad
v rozmezí ±200◦ C je pomˇ ernˇ e pˇ resnˇ e vyjádˇ ren funkcí α(T ) = aT 2 + bT + C kde a = −8.27 × 10−11 , b = 2.21 × 10−8 , c = 1.17 × 10−5 ñ
ñ
koeficienty získány regresí (proložením kˇ rivky) experimentálních dat metodu regrese budeme probírat pozdˇ eji
3/22
PA081: Programování numerických výpoˇ ct˚ u
Montáž ložiska ñ
polomˇ er po zchlazení poˇ cítáme pro ∆T → 0
A. Kˇ renek Montáž ložiska
rx = r0 (1 − ∆T α(T0 ))(1 − ∆T α(T0 − ∆T )) . . . = r0 (1 − ∆T α(T0 ) − ∆T α(T0 − ∆T )
Cirkusové ˇ císlo Obchodní pˇ rípad
+ (∆T )2 α(T0 )α(T0 − ∆T )) . . . ñ
vedlo by na netriviální diferenciální rovnici
ñ
protože α(T ) 1, ∆T 1, lze nelineární ˇ cleny zanedbat
ñ
dostáváme zjednodušení Z T0 ˆ x )) rx = r0 (1 − α(T
kde
ˆ x) = α(T
α(T )dT Tx
a budeme ˇ rešit rovnici ˆ x) = α(T
r0 − rx r0 4/22
PA081: Programování numerických výpoˇ ct˚ u
Montáž ložiska
A. Kˇ renek
ñ
Montáž ložiska
ˇ rešená funkce
Cirkusové ˇ císlo
T0
a 3 b 2 r0 − rx r0 − rx ˆ x )− f (Tx ) = α(T = T + T + cT − r0 3 2 r0 Tx a b r0 − rx = (T03 − Tx3 ) + (T02 − Tx2 ) + c(T0 − Tx ) − 3 2 r0
ñ
Obchodní pˇ rípad
potˇ rebujeme ñ ñ
separaci „správného“ koˇ rene zvolit vhodnou numerickou metodu
5/22
PA081: Programování numerických výpoˇ ct˚ u
Montáž ložiska
A. Kˇ renek
ñ
Montáž ložiska
ˇ rešená funkce
Cirkusové ˇ císlo
T0
a 3 b 2 r0 − rx r0 − rx ˆ x )− f (Tx ) = α(T = T + T + cT − r0 3 2 r0 Tx a b r0 − rx = (T03 − Tx3 ) + (T02 − Tx2 ) + c(T0 − Tx ) − 3 2 r0
ñ
potˇ rebujeme ñ ñ
ñ
Obchodní pˇ rípad
separaci „správného“ koˇ rene zvolit vhodnou numerickou metodu
koeficient u Tx3 je kladný (a bylo záporné) ñ ñ
pro Tx → ∞ také f (Tx ) → ∞ pro Tx → −∞ také f (Tx ) → −∞
5/22
PA081: Programování numerických výpoˇ ct˚ u
Montáž ložiska
A. Kˇ renek
ñ
derivace ˇ rešené funkce
Montáž ložiska
f 0 (Tx ) = −aTx2 − bTx − c ñ
Obchodní pˇ rípad
f (Tx ) má jedno lokální minimum a jedno maximum ñ ñ
v bodech, kde f 0 (Tx ) = 0 numericky stabilní výpoˇ cet? ñ ñ
ñ
ñ
Cirkusové ˇ císlo
druhá pˇ rednáška ;-) √ b = −2.21 × 10−8 , D = 6.60 × 10−8 je OK
pˇ ri nepˇ resném urˇ cení by mohla Newtonova metoda zabloudit, viz dále
vypoˇ ctené extrémy T1 = −265.544 T2 = 532.775
f (T1 ) = 0.000867609 f (T2 ) = −0.00614507
6/22
Montáž ložiska
PA081: Programování numerických výpoˇ ct˚ u A. Kˇ renek
ñ
dva „nesmyslné“ koˇ reny T− < T1 a T+ > T2 ñ ñ ñ ñ
nemají fyzikální význam T+ – ˇ cep se pˇ ri vysokých teplotách zaˇ cne smršt’ovat T− – pˇ ri nízkých teplotách se zaˇ cne znovu roztahovat obojí d˚ usledkem aplikace regresní kˇ rivky mimo rozsah hodnot, ze kterých byla urˇ cena
ñ
skuteˇ cný koˇ ren v intervalu [T1 , T2 ]
ñ
funkce je spojitá, koˇ ren je právˇ e jeden ñ
ñ
Cirkusové ˇ císlo Obchodní pˇ rípad
dokonalá separace
umíme jednoduše poˇ cítat derivaci ñ
ñ
Montáž ložiska
navíc je na (T1 , T2 ) derivace nenulová
Newtonova metoda je ideální
7/22
Montáž ložiska
PA081: Programování numerických výpoˇ ct˚ u A. Kˇ renek Montáž ložiska Cirkusové ˇ císlo Obchodní pˇ rípad
ñ
numerická stabilita výpoˇ ctu f (T ) a f 0 (T ) ñ ñ ñ
ñ
rozsah T ˇ rádovˇ e ve stovkách ˇ cleny polynomu v ˇ rádech 10−5 –10−2 radˇ eji použijeme typ double
potenciální problém Newtonovy metody f 0 (T ) → 0 ñ
nenastane díky známým vlastnostem funkce
8/22
PA081: Programování numerických výpoˇ ct˚ u
Montáž ložiska ñ
Pr˚ ubˇ eh výpoˇ ctu funkcí rtnewt ñ ñ
ñ
A. Kˇ renek
−3
požadovaná absolutní pˇ resnost 10 nemá smysl více, vstupní data (koeficienty a, b, c) jsou ménˇ e pˇ resné a jak budeme pˇ ri stavbˇ e mostu chladit ˇ cep na takhle pˇ resnou teplotu?
T
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ˇ resnosti
ñ
výsledek je -90.056◦ C
Montáž ložiska Cirkusové ˇ císlo Obchodní pˇ rípad
9/22
Cirkusové ˇ císlo
PA081: Programování numerických výpoˇ ct˚ u A. Kˇ renek Montáž ložiska Cirkusové ˇ císlo
Artisté pˇ ripravují nové vystoupení.
Obchodní pˇ rípad
Na zaˇ cátku se jeden z nich zachytí rukama i nohama v kruhové konstrukci o pr˚ umˇ eru 180 cm a další konstrukci vykoulí na plošinu. Je tˇ reba, aby se na konci tohoto pohybu p˚ uvodnˇ e nejvyšší bod kruhové konstrukce dotýkal plošiny a byl ve stejné výšce jako na zaˇ cátku. Jaká je potˇ rebná délka a sklon plošiny?
10/22
PA081: Programování numerických výpoˇ ct˚ u
Cirkusové ˇ císlo
A. Kˇ renek Montáž ložiska Cirkusové ˇ císlo Obchodní pˇ rípad
β α R
Rα β
11/22
PA081: Programování numerických výpoˇ ct˚ u
Cirkusové ˇ císlo
A. Kˇ renek
ñ
Montáž ložiska
oznaˇ címe ñ ñ ñ
Cirkusové ˇ císlo
R – polomˇ er kruhu α – úhel otoˇ cení kruhu (v radiánech) β – úhel naklonˇ ení plošiny
ñ
délka pohybu kruhu po plošinˇ e je Rα
ñ
z naznaˇ ceného ˇ ctyˇ rúhelníku odeˇ cteme ñ ñ
ñ
Obchodní pˇ rípad
α + β = π (další dva úhly jsou pravé) β R tan 2 = Rα
z toho vyplývá ˇ rešená rovnice tan
β 1 = 2 π −β
tedy
f (β) = (π − β) tan
β −1=0 2
12/22
PA081: Programování numerických výpoˇ ct˚ u
Cirkusové ˇ císlo
A. Kˇ renek
ñ
separace koˇ rene ñ ñ ñ ñ ñ
ñ
π
β je ostrý úhel, tj. β ∈ (0, 2 ) ˇ rešení je právˇ e jedno z geometrické podstaty problému formulací rovnic jsme nepˇ ridali falešný koˇ ren dáme si práci s exaktním d˚ ukazem nebo to rovnou zkusíme
Cirkusové ˇ císlo Obchodní pˇ rípad
pro jistotu ovˇ eˇ ríme f (0) = −1
ñ
Montáž ložiska
f (π /2) = 0.5707963
numericky nebezpeˇ cná by byla až oblast β → π ñ ñ
vedla by na souˇ cin typu 0 · ∞ pohybujeme se v bezpeˇ cné vzdálenosti
ñ
budeme pˇ redstírat, že neumíme spoˇ cítat derivaci
ñ
ukážeme metodu seˇ cen, Riddersovu, a Brentovu
13/22
Cirkusové ˇ císlo Metoda seˇ cen, požadovaná absolutní pˇ resnost 10−4
PA081: Programování numerických výpoˇ ct˚ u A. Kˇ renek Montáž ložiska Cirkusové ˇ císlo
n 1 2 3 4 5 6 7 8 9
beta +0.0000000 +1.5707963 +0.7853982 +1.1780972 +0.9817477 +0.8835729 +0.8344855 +0.8099419 +0.8222137
f(beta) -1.0000000 +0.5707963 -0.0240323 +0.3119657 +0.1544612 +0.0679638 +0.0226702 -0.0005027 +0.0111281
15 16
+0.8105171 +0.8104212
+0.0000445 -0.0000467
Obchodní pˇ rípad
...
ñ
poˇ cínaje tˇ retím krokem výpoˇ cet vždy ve stˇ redu intervalu
ñ
pomalá konvergence 14/22
Cirkusové ˇ císlo Riddersova metoda, požadovaná absolutní pˇ resnost 10−4
PA081: Programování numerických výpoˇ ct˚ u A. Kˇ renek Montáž ložiska
n 1 2 3 4 5 6 7 (8)
beta +0.0000000 +1.5707963 +0.7853982 +0.8103685 +1.1905824 +0.8104702 +1.0005263 +0.8104702
f(beta) -1.0000000 +0.5707963 -0.0240323 -0.0000968 +0.3213144 -0.0000001 +0.1704016
Cirkusové ˇ císlo Obchodní pˇ rípad
4. první iteraˇ cní výpoˇ cet 5. p˚ ulení intervalu mezi 2. a 4. 6. další iterace 7. p˚ ulení mezi 5. a 6. 8. výsledek, už je v toleranci 15/22
PA081: Programování numerických výpoˇ ct˚ u
Cirkusové ˇ císlo Brentova metoda, požadovaná absolutní pˇ resnost 10−4
A. Kˇ renek Montáž ložiska Cirkusové ˇ císlo
n 1 2 3 4 5 6 7 (8)
beta +0.0000000 +1.5707963 +1.0000000 +0.7931377 +0.8115179 +0.8104759 +0.8104259 +0.8104759
f(beta) -1.0000000 +0.5707963 +0.1699574 -0.0165737 +0.0009960 +0.0000054 -0.0000422
ñ
poˇ cítá primárním zp˚ usobem, nedošlo na bezpeˇ cné p˚ ulení interval˚ u
ñ
srovnatelnˇ e rychlá konvergence s Riddersem
ñ
pˇ ri vˇ etší požadované pˇ resnosti se liší ± o jeden krok
Obchodní pˇ rípad
16/22
Cirkusové ˇ císlo
PA081: Programování numerických výpoˇ ct˚ u A. Kˇ renek Montáž ložiska Cirkusové ˇ císlo
ñ
všechny metody se v rámci požadované tolerance shodují
ñ
rychlost konvergence dle oˇ cekávání ñ
ñ ñ ñ
Obchodní pˇ rípad
zdvojnásobení požadované pˇ resnosti zdvojnásobí poˇ cet krok˚ u p˚ ulení interval˚ u Riddersova metoda naroste ze 7 na 9 Brentova dokonce jen ze 7 na 8 √ lepší výsledky než teoretická rychlost konvergence 2
17/22
PA081: Programování numerických výpoˇ ct˚ u
Cirkusové ˇ císlo
A. Kˇ renek Montáž ložiska Cirkusové ˇ císlo
ñ
délka pohybu po plošinˇ e je
Obchodní pˇ rípad
Rα = R(π − β) = 2.098 ñ
zbývá spoˇ cítat ˇ cást plošiny od zemˇ e k bodu dotyku tan
π −β R = 2 l
tedy
l=
R = 0.3861 tan((π − β)/2)
ñ
nehrozí významné numerické nepˇ resnosti
ñ
celková délka plošiny je 2.484 m, sklon 46.44◦
18/22
Obchodní pˇ rípad
PA081: Programování numerických výpoˇ ct˚ u A. Kˇ renek Montáž ložiska
Americká garážová firma zahajuje montáž poˇ cítaˇ cu ˚. Kolik jich musí v prvním roce prodat, aby byla zisková?
Cirkusové ˇ císlo Obchodní pˇ rípad
19/22
Obchodní pˇ rípad
PA081: Programování numerických výpoˇ ct˚ u A. Kˇ renek Montáž ložiska
Americká garážová firma zahajuje montáž poˇ cítaˇ cu ˚. Kolik jich musí v prvním roce prodat, aby byla zisková? ñ
Obchodní pˇ rípad
uvažované náklady ñ ñ ñ ñ
ñ
Cirkusové ˇ císlo
fixní jednorázové (poˇ rízení vybavení): $20000 fixní roˇ cní (nájem, web, . . . ): $15000 variabilní (komponenty, hodinová sazba práce, . . . ): $625n semivariabilní (nároˇ cnˇ ejší logistika atd.) $30n1.5
pˇ redpokládané pˇ ríjmy ñ ñ
prodej poˇ cítaˇ cu ˚: $1500n slevy (množstevní, VIP zákazníci, . . . ): −$10n1.5
19/22
PA081: Programování numerických výpoˇ ct˚ u
Obchodní pˇ rípad
A. Kˇ renek Montáž ložiska
ñ
f (n) = T C(n) − T S(n) = 35000 − 875n + 40n ñ
ñ
Obchodní pˇ rípad
=0
prochází bodem (0,35000) osu x protíná v 35000/875 = 40
ˇ clen 40n1.5 zp˚ usobí „prohnutí“ nahoru ñ ñ
ñ
1.5
první dva ˇ cleny – lineární funkce ñ
ñ
Cirkusové ˇ císlo
ˇ rešená rovnice
posune koˇ ren z bodu 40 dál pˇ ridá druhý koˇ ren pro vyšší n
zisku dosahujeme v oblasti mezi tˇ emito koˇ reny
20/22
PA081: Programování numerických výpoˇ ct˚ u
Obchodní pˇ rípad
A. Kˇ renek
ñ
separace 1. koˇ rene ñ ñ
ñ
podle pˇ redchozí úvahy f (40) > 0, skuteˇ cná hodnota 10119 zkusíme f (80) = −6378, vyhovuje
Cirkusové ˇ císlo Obchodní pˇ rípad
separace 2. koˇ rene ñ ñ
pˇ redchozí f (80) = −6378 hrubý odhad zanedbáním absolutního ˇ clene −875n + 40n1.5 = 0
ñ
ñ
Montáž ložiska
tedy
n=
875 40
2 = 478.51
zkusíme f (500) = 44713, vyhovuje
numerická stabilita v dané separaci koˇ ren˚ u ñ ñ ñ
souˇ cty/rozdíly ˇ císel s minimálním ˇ rádovým rozdílem reálnˇ e nás zajímá pˇ resnost na jednotky žádný potenciální problém
21/22
PA081: Programování numerických výpoˇ ct˚ u
Obchodní pˇ rípad
A. Kˇ renek
ñ
vypoˇ ctené koˇ reny 62.7 a 384.0
Montáž ložiska
ñ
výsledky v požadované pˇ resnosti shodné všemi metodami
Cirkusové ˇ císlo
ñ
chování jednotlivých metod (poˇ cty iterací) pˇ resnost koˇ ren p˚ ulení Ridders Brent
10−1 1 2 11 15 5 9 5 9
10−4 1 2 21 25 7 11 7 10
ñ
umˇ elá pˇ resnost 10−4 – zdvojnásobení poˇ ctu platných ˇ císlic
ñ
dokládá oˇ cekávanou rychlost konvergence ñ ñ
Obchodní pˇ rípad
lineární pro p˚ ulení interval˚ u √ 2 pro ostatní metody
22/22