Aalto-yliopisto, Matematiikan ja Systeemianalyysin laitos Apiola/Kuortti

Mat-1.C Matemaattiset ohjelmistot

Harjoitus 2 18-21.3.2013

• Kaikki 2-tehtävät palautetaan ke 27.3.2013 klo 24.00 mennessä Petrille osoitteeseen: petri.leskinen@aalto.fi. Laita liitteeksi m-tiedostot ja Maple-ws:t (.mw), poista jälkimmäisistä mielellään tulostukset Edit-valikon “Remove output”-valinnalla.

• Merkitse selvästi, mitkä ovat varsinaisia DOKU-tehtäviä. Hyvä, jos teet mailiviestiin luettelon liitteinä lähettämiesi tehtävien numeroista, ja merkitset siihen D:t.

• Ma 25.3. opettajat paikalla (niin kauan kuin oppilaita riittää), läsnäoloa ei tarkisteta (syy: välikokeet).

Maple-työt tehdään työarkille, lataa Mapleen pohja.mw. Käy läpi kokeiluja ohjatusti tehden:

maplepohjanew.mw, mapleperusteita.mw

Hyödyllisiä vihjeitä:

• Kursoria ei tarvitse siirtää rivin loppuun ennen Enter-käskyä!

• Nuolinäppäimillä voi siirtyä yläindeksistä pois; samoin murtolausekkeissa.

• Pikanäppäimiä:

Ctrl + Delete poistaa käsky- tai tulosrivin

Ctrl + t siirtyy tekstitilaan

F5 siirtyy tekstitilassa kaavankirjoitustilaan ja takaisin Ctrl + k tekee uuden käskyrivin kursorin yläpuolelle Ctrl + j tekee uuden käskyrivin kursorin alapuolelle Ctrl + l (l = label) liittää viittauksen aikaisemman tuloksen numeroon

Lisää Maple-ohjeita lopussa.

1. mplDi0001.tex Maple Määritä funktion p(x) = x^{3} − 4x^{2} + 4x − 1 nollakohdat ja lokaali max/min. Piirrä funktion ja derivaatan kuvaajat.

Vihje: solve,evalf,diff,plot. Yksinkertaisinta ehkä käsitellä lausekkeena, mutta saat kokeilla myös funktiotapaa.

2. mplDi010.tex (Harj. 1:n Matlab-tehtävä nyt Maplella) Määritä funktion f(x) = arcsin(2x\sqrt{1 − x^{2}} ) suurin ja pienin arvo välillä [−1, 1].

Käytä symboliohjelmissa perinteistä “diffistekniikkaa” kuvan kanssa, Matlab:ssa raakaa “numeronmurskausta” tyyliin: linspace, plot, zoom, uusi linspace kapeammalla välillä, find, ...

Vihje: arcsin on Mathematicassa ArcSin, Maplessa arcsin ja Matlabissa asin. Symbolilaskentaohjelma saattaa johtaa oikeaan tulokseen puutteellisin perustein, jos tarkkoja ollaan.

3. mplDi011.tex Laske sen alueen pinta-ala, jota rajoittavat käyrät y^{2} = x ja x − y = 3. Kuva, tottakai!

Vihje: Mieti, kumpi on helpompaa: integrointi x- vai y-suunnassa.

4. mplDi017.tex

a) Osoita, että funktio arctan \frac{y}{x} toteuttaa Laplacen osittaisdifferentiaaliyhtälön

\frac{∂^{2}u}{∂x^{2}} + \frac{∂^{2}u}{∂y^{2}} = 0.

(Tällaisia funktioita sanotaan harmonisiksi funktioiksi.)

b) Oletetaan, että funktioilla u(x, y) ja v(x, y) on jatkuvat toiset osittaisderivaatat ja ne toteuttavat ns. Cauchy-Riemannin yhtälöt:

\frac{∂u}{∂x} = \frac{∂v}{∂y}, \frac{∂v}{∂x} = −\frac{∂u}{∂y}

Osoita, että u ja v ovat harmonisia.

c) Olkoon f(x, y) = x^{3}y^{2} + x^{4} siny + cos(xy). Laske osittaisderivaatat f_{xxy}, f_{xyx}, f_{yxx} ja totea, että ne ovat samat.

5. mplDi005.tex

Laske integraali

∫_{0}^{2π} \frac{cosx}{13 − 12 cos2x} dx

a) symbolisesti, b) numeerisesti. Piirrä integroitavan funktion kuvaaja. Mikä itse asiassa on integraalin arvo?

Suorita b)-kohta sekä Maplella että Matlab:lla (Toki voit kokeilla Matlab:lla myös a)-kohtaa tyyliin syms x ja sitten Maple-komento int.)

Vihje: Mathematica: Symbolinen integrointi tapahtuu funktiolla Integrate, numeerinen funktiolla NIntegrate. Jälkimmäisessä sovelletaan suoraan jotakin numeerisen integroinnin menetelmää, jonka valintaan myös käyttäjä voi vaikuttaa. Ks. dokumentaatiota, erityisesti Implementation Notes.

Maple: Symbolinen integrointi tapahtuu funktiolla int, numeerinen funktiolla int(...,type=numeric) tai evalf(Int(...)). Numeerisessa sovelletaan suoraan jotakin numeerisen integroinnin menetelmää, jonka valintaan myös käyttäjä voi vaikuttaa.

Esim: evalf(Int(f, x = 0 .. 2, digits = 20, method = _Dexp))

Matlab: Integrandi määritellään funktioksi (helpoimmin funktiokahvaksi “function handle”). Sitten quad-alkuiset Matlab-funktiot.

Viitteet:

http://math.tkk.fi/~apiola/matlab/opas/lyhyt/m-files.html (Matlab:n funktiokahva,

function handle)

6. mplDi0011.tex Kuulantyönnön tulos riippuu kuulan alkunopeudesta v, lähtökorkeudesta h ja työnnön suuntakulmasta x seuraavan lausekkeen mukaisesti:

f := x → \frac{v cos(x) (v sin(x) + \sqrt{v^{2} sin(x)^{2} + 2hg} )}{g} ,

missä x ∈ [−π/2, π/2]. Käytetään SI-järjestelmän yksiköitä ja oletetaan, että h = 2, v = 14 ja g = 9.81. Määritä työnnön optimaalinen suuntakulma ja maksimitulos.

Kannattanee edetä seuraavien vaiheiden mukaan:

• Määrittele f funktiona; älä sijoita lukuarvoja tässä vaiheessa, niin voit tarkistaa, että lauseke on oikein.

• Sijoita lukuarvot h, v, g.

• Piirrä funktion f kuvaaja välillä −π/2 ≤ x ≤ π/2 ja tarkista, että se näyttää järkevältä. (Yleinen virhe: kertomerkkejä puuttuu!)

• Ratkaise maksimi kokeilemalla molempia tapoja: suoraan maximize TAI muodosta yhtälö f′(x) = 0, ratkaise numeerisesti fsolve-käskyllä, laske maksimi.

• Muuta saatu kulma asteiksi ja mieti, onko tulos järkevä.

7. mplD0009.tex Putoavan kappaleen nopeus v = v(t) toteuttaa differentiaaliyhtälön mv′(t) = mg − kv(t)^{2}, jos positiivinen suunta on alaspäin ja ilmanvastus on verrannollinen nopeuden neliöön kertoimella k > 0.

a) Ratkaise differentiaaliyhtälö alkuehdolla v(0) = 0.

b) Mikä on rajanopeus lim_{t→∞} v(t)? Vihje: Ohjelma ei osaa laskea raja-arvoa, koska se ei tiedä vakioiden etumerkkiä. Lisää käsky assume(m>0 and k>0 and g>0) ja kokeile uudelleen sen jälkeen.

8. mplDi0006.tex a. Olkoon f(x, y) = \frac{x−y}{x+y} .

Määritä pinnan z = f(x, y) tangenttitaso pisteessä (2, −1). Piirrä pinta ja tangenttitaso ja pyörittele ja zoomaa.

b. Sama pinnalle z = arctan \frac{y}{x} pisteessä (2, 2, π/4).

9. mplDi0002.tex Osittaisderivoituvan funktion f : ℝ^{2} → ℝ gradientti ∇f määritellään näin ∇f(x, y) = f_{x}(x, y)i + f_{y}(x, y)j .

Olkoon f(x, y) = |xy| .

(a) Piirrä tasa-arvokäyrät(korkeuskäyrät) f(x, y) = k, k = 1, 2, 3.

(b) Piirrä f:n gradienttivektoreita ∇f(x, y) tasa-arvokäyrien pisteisiin. Kun käytät samaa skaalaa akseleilla (scaling=constrained), pitäisi kuvasta näkyä, miten gradienttivektorin ja korkeuskäyrän suunnat suhtautuvat toisiinsa.

Vihje: Aloita työarkki näin:

> restart:
> with(plots): with(plottools):
> nuoli:=(alkup, loppup,vari)->arrow(alkup,loppup,0.01,0.05,0.02,color=vari);
> korkeuskayra:=k->implicitplot(abs(x*y)=k,x=-2..2,y=-2..2);
> #  Maariteltiin grafiikka-arvoinen funktio, usein tosi katevaa!
> kkparvi:=display(seq(korkeuskayra(k),k=1..3);
>

Kts. lisää: mplDi0002Apu.mw

10. mplDi0004.tex Maaston korkeus (merenpinnasta mitattuna) karttakoordinaattien funktiona olkoon

h(x, y) = −x^{2} + 4xy − 8y^{2} + 300.

Positiivinen x-akseli osoittaa itään ja positiivinen y-akseli pohjoiseen.

a. Kulkuri K ottaa pisteestä (1, 2, h(1, 2)) lähtöaskeleen kaakkoon. Nouseeko hän vai laskeutuuko? Tämä on käsinlaskutehtävä, mutta tee Maplella. Havainnollista Maplepiirroksin: Pintapiirros: plot3d, korkeuskäyrät: contourplot tai implicitplot. Leikkauskäyrä kaakko-luode-suuntaisen pystytason kanssa.

b. Muodosta funktion h(x, y) gradienttifunktio (gradienttikenttä). Piirrä gradienttikenttä plots-pakkauksen funktiolla fieldplot. Yhdistä korkeuskäyräpiirros tämän kanssa display-funktion avulla.

Vihje: Gradienttikentän voi laskea (tietysti käsin) tai derivoimalla Maplen diff:llä tai linalg-pakkauksen funktiolla grad. Ei ole pahitteeksi, jos kokeilet kaikkia tapoja.

11. DOKU mplV0001.tex Joudut tekemään vastuunalaisen päätöksen mitoista valmistettaessa laatikkoa. Pohjamateriaali on kaksi kertaa niin kallista pinta-alayksikköä kohti kuin sivu- tai kansimateriaali. Millä mitoilla saat V-tilavuuksisen laatikon materiaalikustannukset minimoiduksi? Perustele, että ratkaisusi on globaali minimi joukossa {(x, y)|x > 0, y > 0}. (Toisia derivaattoja ei välttmättä tarvita.) (Käytä Maplea sydämesi kyyllydestä, vaikka käsinlaskullakin selviäisit.)

12. mplDi0005.tex Määritä lieriöiden

x^{2} + y^{2} = 2, y^{2} + z^{2} = 2

leikkauskäyrän pisteen (1, −1, 1) kautta kulkevan tangentin yhtälö.

Lieriöpinnan piirtäminen sujuu hyvin plot3d:llä. Kannttaa ajatella lieriö (kahdesta parametristä riippuvana) parametrimuotoisena pintana. Ensimmäisen lieriön luonnollinen parametriesitys on x = \sqrt{2}cost, y = \sqrt{2}sint, z = z. Tässä siis t ja z ovat parametreja.

Edellinen voisi näyttää tältä:

plot3d([sqrt(2)*cos(t),sqrt(2)*sin(t),z],t=0..2*Pi,z=c..d);

Jälkimmäinen vastaavasti. Kuvat yhdistetään:

display(kuva1,kuva2);

Huom! plot3d on monipuolinen funktio, sille voi antaa pinnan muodossa f(x, y), mutta myös parametrimuodossa yllä kaavailtuun tapaan.

13. DOKU mplV0009.tex Määritä funktion

f(x, y) = 1/x + 1/y + sin(x^{2}y^{2})

suurin ja pienin arvo joukossa [1, 2] × [1, 2].

Maplen lisäksi kannattaa kokeilla Matlab:lla meshgrid, max/min, find ...-tekniikkaa. Toki ihan optimointiin räätälöityjä funktioitakin kummallakin ohjelmalla. Mutta ensisijaisesti ihan perustekniikoita, please!.

Vihje: Vaatimattomasta ulkoasustaan huolimatta voi olla hiukan työläs.

14.

DOKU Maple,Matlab (H2T10)

a) Ratkaise alkuarvotehtävä

y′− y = cosx, y(0) = 1

analyyttisesti Maplella ja numeerisesti Matlabilla. Piirrä ratkaisukäyrä. b) Anna alkuarvoksi symboli c ja piirrä ratkaisukäyräparvi sopivalla välillä, kun c =

−0.9, −0.8, …, 0.

Miltä parvi näyttää suurilla x :n arvoilla. Tässä pitäisi erottua kolmenlaista käytöstä.

Vihje: Maple: dsolve, Matlab: ode45

Avainsanat: Differentiaaliyhtälö, alkuarvotehtävä, analyyttinen ratkaisu, numeerinen ratkaisu.

15. DOKU Ratkaise differentiaaliyhtälöryhmä

\{ \begin{matrix}\frac{dx}{dt} = −σx + ρy \\ \frac{dy}{dt} = σx − y − xz \\ \frac{dz}{dt} = −βz + xy \end{matrix}

numeerisesti välillä [0, 20], kun σ = 10, ρ = 28 ja β = 8/3. Piirrä ratkaisukäyrät samaan kuvaan, ja piirrä käyrät x(t) ja z(t) parametrisesti. Tämän jälkeen piirrä 3-ulotteinen parametrisoitu käyrä kaikista koordinaateista.

Onko ratkaisu rajoitettu? Suppeneeko se kohti jotain arvoa?

Kokeile muuttaa alkuarvoja, sekä parametrien arvoja. Vallitsevan teorian mukaan systeemi on kaoottinen dynaaminen systeemi, jonka käyttäytyminen voi muuttua merkitsevästi jo pienistä muutoksista lähtötilanteessa; itse asiassa termi perhosvaikutus keksittiin kuvaamaan juuri tämän systeemin käytöstä.

Vihje: Kolmiulotteinen parametrisoitu käyrä (tai pistejoukko) piirretään MATLABissa funktiolla plot3.

16. mplV000.tex

Ohjeita

Kerätään ohjeita näiden tehtävien aihepiiriin liittyen. “Tehtävä”-linkistä saat L^{A}TEX-koodin, josta sopivan osan voit haluamallasi tavalla muokaten liittää tehtäväpaperiisi.

Neliömuotojen definiittisyys

Määr: Neliömuoto q(x) = x^{T}Ax (A on symmetrinen matriisi) on

1. positiivisesti definiitti, jos q(x) > 0 ∀x ≠ 0,

2. negatiivisesti definiitti, jos q(x) < 0 ∀x ≠ 0,

3. positiivisesti semidefiniitti, jos q(x) ≥ 0 ∀x ∈ ℝ^{n} ja ∃y ≠ 0, jolla q(y) = 0, 4. negatiivisesti semidefiniitti, jos q(x) ≤ 0 ∀x ∈ ℝ^{n} ja ∃y ≠ 0, jolla q(y) = 0,

5. indefiniitti, jos ∃x, y siten, että q(x) > 0 ja q(y) < 0.

Samoja definiittisyyskäsitteitä käytetään myös symmetrisestä matriisista A.

Suunnattu derivaatta ja gradientti

• Suunnattu derivaatta pisteesä p_{0} vektorin v suunassa saadaan lasketuksi pisteessä p_{0} lasketun gradientin ja suuntayksikkövektorin sisätulona.

• Siispä funktio kasvaa nopeimmin gradientin suuntaan ja sen kasvu on 0 gradienttia vastaan kohtisuoraan suuntaan.

• Suunta, johon funktion kasvu on 0 on tasa-arvokäyrän (tai -pinnan) tangentin (tangenttitason) suuntainen, joten gradientti on normaalin suuntainen.

Pinnan normaali ja tangenttitaso

Jos pinnan yhtälö esitetään muodossa F (x, y, z) = 0, saadaan edellisen perusteella pinnan tangenttitason yhtälö pisteessä p_{0} näin:

∇F (p_{0})˙(p − p_{0}) = 0

Jos pinta on annettu muodossa z = f(x, y), saadaan siten normaalin suunta funktion F (x, y, z) = f(x, y) − z gradienttina.

Tästä seuraa, että pisteeseen p_{0} asetetun tangenttitason yhtälö voidaan kirjoittaa muotoon

z − z_{0} = f_{1}(p_{0})(x − x_{0}) + f_{2}(p_{0})(y − y_{0}).

(f_{1}ja f_{2} tarkoittavat osittaisderivaattoja.)

Kahden pinnan leikkauskäyrän tangentti

Leikkauskäyrän tangentti on kohtisuorassa molempien pintojen normaalia vastaan (eikö vain!). Siten leikkauskäyrän tangentin suuntainen vektori saadaan pinnan normaalivektorien ristitulona t = n_{1} × n_{2}.

Kriittiset pisteet, ääriarvot

Kriittinen piste (KRP) p: ∇f(p) = 0. Kriittisen pisteen laatu selviää (jos selviää) Hessen matriisin H_{f} (p) definiittisyydestä. Symmetrisen matriisin definiittisyyskäytös selvitetään ominaisarvojen avulla. Jos matriisi on 2⋊2, voidaan käyttää determinanttia (kts. tehtävä mplV006a). Isommillekin matriiseille on determinanttiehtoja, mutta ne on hankala muistaa ja käyttää, jääkööt muistoksi “determinanttien kulta-ajoilta”.

Hessen matriisi saadaan jälleen yksinkertaisimmin tähän tapaan:

 > with(linalg)
 > H:=Matrix(hessian(f,[x,y]));
 > H00:=subs(x=0,y=0,H)  % Arvo pisteessä (0,0)

(Tietysti on helppo kirjoittaa ihan oma “hessi”.)

Vektorikenttä ja gradientti

with(linalg): with(plots):
fieldplot(grad(f(x,y),[x,y]),x=a..b, y=c..d,arrows=slim,color=x);
# a:lla, b:lla jne. oltava tietysti numeeriset arvot.

Uusissa Maplen versioissa on kirjastopakkaus VectorCalculus ja siellä funktio Gradient lukuisine valitsimineen. Kts. helppi. Vanhan linalg-kirjaston kunnon grad on perustarpeisiin ehkä yksinkertaisin ja helppokäyttöisin.

Oma pikku funktio on usein selkein, se voidaan määritellä ongelmakohtaisesti esim. toimimaan vain 2d-tilanteessa. Tällainen gradienttifunktio voitaisin kaikessa yksinkertaisuudessaan

määritellä näin: gradi2:=(f,x,y)->[diff(f,x),diff(f,y)]

Pintapiirroksen “valaiseminen” esim. avaruuskäyrillä

Usein pintapiirrosta voidaan täsmentää ja tarkentaa ja ymmärtää paremmin, kun piirretään sopivia avaruuskäyriä pinnalle spacecurve:lla.

Alla on esimerkki, jossa piirretään napasädettä pitkin kulkevan pystytason ja annetun pinnan leikkauskäyrä sekä käyrän projektio xy-tasossa. Varsin käyttökelpoinen tapa monessa yhteydessä. Tätä voi modifioida tarpeen mukaan.

> with(plots):
> f:=(x,y)->4-x^2-y^2;
> x:=r*cos(Theta):y:=r*sin(Theta): Theta:=Pi/4:
> pystyleikkaus:=spacecurve({[x,y,f(x,y)],[x,y,0]},r=0..2,thickness=3,
color=blue,axes=BOX)
> x:='x':y:='y':  # On hyv\"a muistaa vapauttaa.
> pinta:=plot3d(...):   # Muista t\"ass\"a tavassa lopettaa kaksoispisteeseen.
> display([pinta,pystyleikkaus],style=patchcontour);
> display(pystyleikkaus);  # Katsotaan pelkk\"a\"a leikkausk\"ayr\"a\"a.