> restart; # Egyenletek, egyenletrendszerek megoldása # # Sokszor van szükség egyenletek, egyenletrendszerek megoldására, ide # beleértve # lineáris egyenletrendszerek, illetve differenciálegyenletek megoldását # is. # Gyakran a megoldást nem tudjuk előállítani pontosan, így meg kell # elégedni # ilyenkor valamilyen közelítő megoldással. Ezen a munkalapon bemutatunk # néhány egyszerű lehetőséget a fenti problémák kezelésére. # # 1. rész: Egyenletek megoldása # # Az egyenlet két kifejezés logikai egyenlővé tételéből áll. # Megadásánál az = jelet használjuk a két oldal között. Vigyázat! Ne # keverjük ezt össze az értékadás := jelével! # Az egyenlet bal oldalára az lhs, jobb oldalára az rhs paranccsal # hivatkozhatunk. # # Pl. > egy1:=x^3+x^2-2*x+4=24/x; 3 2 egy1 := x + x - 2 x + 4 = 24 1/x > lhs(egy1); 3 2 x + x - 2 x + 4 > rhs(egy1); 24 1/x # Az egyenletből el is hagyható az = jel és a jobb oldal; ilyenkor a # Maple automatikusan hozzáérti a 0-át jobb oldalként. # Pl. > egy2:=exp(2*x^2)-4*cos(x); 2 egy2 := exp(2 x ) - 4 cos(x) # Ha már van egy egyenletünk, akkor ezt megoldatni a solve paranccsal # lehet. # Szintaxisa: solve(egyenlet, változó); > solve(egy1, x); 2, -3, 2 I, -2 I # Vagy egy sorban: > solve(x^3+x^2-2*x+4=24/x, x); 2, -3, 2 I, -2 I # Általában az egyenlet gyökei nem ilyen szépek. Ilyenkor felbukkan a # RootOf # parancs, jelezve, hogy az utána zárójelben következő kifejezés gyökeit # kell a képletbe # helyettesíteni. Ezek mindegyikének kiértékelése pl. az allvalues # utasítással lehet. # Pl. > RootOf(exp(4*x)-3*exp(2*x)=1/2, x); RootOf(2 exp(4 _Z) - 6 exp(2 _Z) - 1) > allvalues(%); 1/2 ln(3/2 + 1/2 sqrt(11)), 1/2 ln(3/2 - 1/2 sqrt(11)) # Néha a Maple nem mond semmit, azt jelezve, hogy nem talált gyököt. # Pl. a fenti 2. egyenletnél: > egy2; #csak emlékeztetőül kiiratva 2 exp(2 x ) - 4 cos(x) > solve(egy2, x); # Tényleg nincs gyök? > plot(egy2, x=-1..1); # Dehogy nincs! Van, de a Maple zárt alakban ezt nem tudja megadni. # A közelítő megoldásra az fsolve parancs szolgál, harmadik # argumentumként megadva azt az intervallumot, ahol a gyököt keressük: > fsolve(egy2, x, 0.6..0.8); .7368730662 > fsolve(egy2, x, -0.8..-0.6); -.7368730662 # Rábízhatjuk a Maple-re, hogy keresse meg saját maga a megfelelő # intervalumot, de rendszerint így nem kapjuk meg az összes gyököt: > fsolve(egy2, x); .7368730662 # Ebből az jönne ki, hogy csak egy gyök van, ami hibás. Tehát rajzoljunk # grafikont és gondolkodjunk, mielőtt az fsolve-ot használjuk! # # Két megjegyzés. # 1. Gyakran a Maple segédfüggvényekkel adja meg a gyököt. Ha ez nem # tetszik, # akkor közelítsük az eredményt evalf-fal. # Pl. > egy3:=exp(x)=2-2*x; egy3 := exp(x) = 2 - 2 x > solve(egy3, x); -LambertW(1/2 exp(1)) + 1 > evalf("); # Elenőrzés: > plot(lhs(egy3)-rhs(egy3), x=-10..4); > # #A rajzból az látszik, hogy az eredmény elfogadható. # 2. megjegyzés. Ahogy korábban volt róla szó, az egyenlet tartalmazhat # paramétereket is. # Pl. > solve(x^2+a=3*x, x); 3/2 + 1/2 sqrt(9 - 4 a), 3/2 - 1/2 sqrt(9 - 4 a) > # Gyakorló feladatok ehhez a részhez. # 1. feladat. # Legyen f az f(x)=x^2 tg(x+1) hozzárendeléssel és a D_f=[-2,2] # értelmezési tartománnyal megadva. # # Állapítsd meg f lokális szélsőértékhelyeit! (Oldd meg: f'(x)=0, # stb.) # # 2. feladat. # Nézd meg, mit ad a Maple a sin(x)=1/2 egyenletre! Találj meg legalább # három gyököt! # # # 2. rész: Egyenletrendszerek megoldása # # Ugyanúgy, mint az egyenleteké, csak most az egyenleteket, illetve az # ismeretleneket # kapcsos zárójelek között fel kell sorolni. # # Néhány példa: # # 1. példa. > solve( {x^2/y=-4, x*y-1=y}, {x,y} ); 3 2 2 3 2 {y = - 1/4 RootOf(_Z + 4 - _Z ) , x = RootOf(_Z + 4 - _Z )} > evalf(allvalues(")); > {x = -1.314596213, y = -.4320408008}, { x = 1.157298106 - 1.305151527 I, y = .09102040058 + .7552246950 I}, { x = 1.157298106 + 1.305151527 I, y = .09102040058 - .7552246950 I} # 2. példa. # Hol van az f(x,y)=(x+2)/(2+x^2+y^2-2*y) fv-nek lokális # szélsőértékhelye? > f:=(x+2)/(2+x^2+y^2-2*y); x + 2 f := ----------------- 2 2 2 + x + y - 2 y > solve( {diff(f,x), diff(f,y)} , {x,y}); 2 {y = 1, x = RootOf(-1 + _Z + 4 _Z)} > evalf(allvalues(%), 5); {x = .2361, y = 1.}, {y = 1., x = -4.2361} # Ellenőrizzük a pontokat pl. így: > plot3d(f, x=-4.5..-4, y=0..2, > axes=boxed,style=patchcontour,shading=zhue, title=`második gyök`); # Látjuk, hogy lok. min. # > plot3d(f, x=0..1, y=0..2, axes=boxed, > style=patchcontour,shading=zhue,title=`első gyök`); # Látjuk, hogy lok. max. # # # 3. rész: Lineáris egyenletrendszerek megoldása # # A megfelelő lineáris algebrai ("mátrixos") Maple-utasításokhoz be kell # olvastatni a Maple linalg csomagját: # > with(linalg): # Mátrixokat pl. sorfolytonos alakban lehet megadni, a szögletes # zárójelek segítségével. # Pl. > A := matrix( [[0,2,3],[-1,2,4]] ); [ 0 2 3] A := [ ] [-1 2 4] > b := matrix( [ [10], [20], [30] ] ); [10] [ ] b := [20] [ ] [30] # Rövidebben is be lehet írni b-t a vector utasítással: > b_vek := vector( [10,20,30] ); b_vek := [10, 20, 30] # Szemre b és b_uj különbözik, de valójában egyenlők. # # Az elemekre való hivatkozás módja: pl. > b_vek[2]; 20 > b[2,1]; 20 > A[2,3]; 4 > # A mátrixműveletek elvégzéséhez jelöljük ki a kívánt műveletet és # rakjuk be az evalm utasítás argumentumába (evalm=evaluate matrix). # A műveleti jelekhez csak annyit, hogy a mátrix-szorzás jele &*. # # Pl. > A := matrix( [[1,2,3], [-3,2,9], [5,4,5]] ); [ 1 2 3] [ ] A := [-3 2 9] [ ] [ 5 4 5] > B := matrix( [[1,2,3], [15,24,35]] ); [ 1 2 3] B := [ ] [15 24 35] > C := matrix( [[1,2], [15,35], [-1,4]] ); [ 1 2] [ ] C := [15 35] [ ] [-1 4] > evalm(2*A-3*A^2); [ -28 -50 -102] [ ] [-114 -98 -144] [ ] [ -44 -106 -218] > evalm(A &* B); Error, (in linalg[multiply]) non matching dimensions for vector/matrix product > evalm(B &* A); [ 10 18 36] [ ] [118 218 436] > evalm(inverse(A) &* C); [ 33 ] [-2/7 -- ] [ 14 ] [ ] [ -151] [-18/7 ----] [ 14 ] [ ] [ 99 ] [15/7 -- ] [ 14 ] # Persze az elemek paramétereket is tartalmazhatnak. Pl. > A := matrix( [[1,2,3], [-3,2,x], [1,4,5]] ); [ 1 2 3] [ ] A := [-3 2 x] [ ] [ 1 4 5] > inverse(A); [ -5 + 2 x 1 x - 3 ] [ -------- - ----- - ----- ] [ 1 + x 1 + x 1 + x ] [ ] [ x + 15 1 x + 9] [- 1/2 ------ - ----- 1/2 -----] [ 1 + x 1 + x 1 + x] [ ] [ 1 1 1 ] [ 7 ----- ----- -4 ----- ] [ 1 + x 1 + x 1 + x ] > det(A); # A determinánsa -2 - 2 x # Egy példa egyenletrendszerek megoldására: # # A megfelelő mátrixok bevezetésével oldd meg az alábbi # egyenletrendszert! # 2y+3z=3 , -x+2y+4z=-2, 2x+y=1 > A := matrix( [[0,2,3],[-1,2,4],[2,1,0]] ); [ 0 2 3] [ ] A := [-1 2 4] [ ] [ 2 1 0] > b := vector( [3,-2,1]); b := [3, -2, 1] > evalm( inverse(A) &* b ); [-16, 33, -21] # Tehát a megoldás: x=-16, y=33, z=-21. # # Általánosabb módszer lineáris egyenletrendszerek megoldására a # linsolve utasítással történik. Ez akkor is működik, ha pl. az # egyenletrendszernek végtelen sok megoldása van (az ezek leírásához # szükséges paraméterek jelzéséhez van szükség linsolve negyedik # argumentumára). # # Az előző egyenletrendszer megoldása: > linsolve(A, b, 'r', t); [-16, 33, -21] # Látjuk, hogy t-re nem volt szükség, mert csak egy megoldás van. Az "r" # változó # az egyenletrendszer rangját adja meg (azt a számot, hogy lényegében # hány különböző # ("hány független") egyenlet van az egyenletrendszerben. # # Egy példa, amikor végtelen sok gyök van: > A := matrix([ [1,2,3], [-2,3,1] ]); [ 1 2 3] A := [ ] [-2 3 1] > b := vector([10,20]); b := [10, 20] > linsolve(A,b,'r',t); [-10/7 - t[1], 40/7 - t[1], t[1]] # Tehát végtelen sok megoldás van: t1 minden választására van egy-egy. > restart; # # 4.rész: Differenciálegyenletek megoldása # # Erre a feladatra - fontosságának megfelelően - a Maple számos # lehetőséget kínál. # Mi most a legegyszerűbbel ismerkedünk meg. # # A megfelelő utasítás: dsolve. Használata hasonló a solve-éhoz: # dsolve(mit, ismeretlen). # Fontos, hogy az ismeretlen függvényt "y(x)-es alakban" írjuk be. # # Egy példa: > de1 := diff(y(x), x)=cos(x)+y(x); d de1 := -- y(x) = cos(x) + y(x) dx > dsolve(de1, y(x)); y(x) = - 1/2 cos(x) + 1/2 sin(x) + exp(x) _C1 # Látjuk, hogy az integrálási állandót "_C1"-gyel jelöli (ha több van, # akkor "_C2", stb.). # # Figyelem! Az eredmény formálisan egy logikai egyenlőség, ennek jobb # oldala szolgáltatja a differenciálegyenlet megoldásának képletét! # # Így a megoldás képlete: > rhs("); > Warning, incomplete string; use " to end the string (note that the ditto operator is now % instead of ") # Kezdeti feltételt is adhatunk meg: > kezdfelt1 := y(Pi/2)=-1; # Ekkor a megoldás: > dsolve({de1, kezdfelt1}, y(x)); exp(x) y(x) = - 1/2 cos(x) + 1/2 sin(x) - 3/2 --------------------------- cosh(1/2 Pi) + sinh(1/2 Pi) # Ennek kirajzolása pl. a [2,3] intervallumon és helyettesítési értéke # x=4-ben 3 jegyre kerekítve: exp(x) y(x) = - 1/2 cos(x) + 1/2 sin(x) - 3/2 --------------------------- cosh(1/2 Pi) + sinh(1/2 Pi) > megold1:=dsolve({de1, kezdfelt1}, y(x)); > y1:=rhs(megold1); exp(x) y1 := - 1/2 cos(x) + 1/2 sin(x) - 3/2 --------------------------- cosh(1/2 Pi) + sinh(1/2 Pi) > plot(y1, x=2..3); > evalf(subs(x=4, y1), 3); kezdfelt1 := y(1/2 Pi) = -1 exp(x) y(x) = - 1/2 cos(x) + 1/2 sin(x) - 3/2 --------------------------- cosh(1/2 Pi) + sinh(1/2 Pi) megold1 := y(x) = exp(x) - 1/2 cos(x) + 1/2 sin(x) - 3/2 --------------------------- cosh(1/2 Pi) + sinh(1/2 Pi) exp(x) y1 := - 1/2 cos(x) + 1/2 sin(x) - 3/2 --------------------------- cosh(1/2 Pi) + sinh(1/2 Pi) -17.2 # # Másik példa: forrás nélküli LRC-kör: > de_LRC:=L*diff(q(t), t,t) + R*diff(q(t),t) + q(t)/C =0; / 2 \ |d | /d \ q(t) de_LRC := L |--- q(t)| + R |-- q(t)| + ---- = 0 | 2 | \dt / C \dt / > dsolve(de_LRC, q(t)); # # Kezdeti feltételt is adhatunk meg: > kezdfelt_LRC:=q(0)=1, D(q)(0)=0; > q_LRC:=rhs( dsolve( {de_LRC,kezdfelt_LRC}, q(t)) ); # Konkrét L,R,C értékekkel: > L:=1; R:=0; C:=1; > #L:='L': R:='R': C:='C': > q_LRC; # Látjuk, hogy az eredményben sqrt(-4) is szerepel. Persze minden valós # végülis, pl. > plot(q_LRC, t=0..10); # Az Euler-formula felhasználásával a komplex alakot átalakíthatjuk # valóssá: > evalc(q_LRC); # # Közelítő megoldásra is szükség lehet, hiszen a megoldás nem mindig # írható fel zárt alakban. # Pl. > de2 := diff(y(x),x)=cos(x^2+y(x)^2); 2 (C R - sqrt(-C (-R C + 4 L))) t q(t) = _C1 exp(- 1/2 --------------------------------) C L 2 (C R + sqrt(-C (-R C + 4 L))) t + _C2 exp(- 1/2 --------------------------------) C L kezdfelt_LRC := q(0) = 1, D(q)(0) = 0 / | q_LRC := 1/2 | \ 2 \ (C R - sqrt(-C (-R C + 4 L))) t | (C R + sqrt(%1)) exp(- 1/2 --------------------------------)| C L / / | /sqrt(%1) - 1/2 | \ 2 \ (C R + sqrt(-C (-R C + 4 L))) t | (C R - sqrt(%1)) exp(- 1/2 --------------------------------)| C L / /sqrt(%1) 2 2 %1 := C R - 4 C L L := 1 R := 0 C := 1 1/2 exp(1/2 sqrt(-4) t) + 1/2 exp(- 1/2 sqrt(-4) t) cos(1/2 sqrt(4) t) d 2 2 de2 := -- y(x) = cos(x + y(x) ) dx > dsolve({de2,y(0)=0}, y(x)); # Nem kaptunk semmit. Ezért közelítő, más néven numerikus megoldáshoz # adjuk be dsolve harmadik argumentumának azt, hogy numeric. Így egy # eljárást # kapunk, amellyel minden x-re ki tudjuk számolni a közelítő megoldást # x-ben. # Az előző példa: > dsolve({de2,y(0)=0}, y(x), numeric); proc(rkf45_x) ... end # Nevezzük ezt el valaminek a könnyebb használhatóság miatt: > nummego2:=dsolve({de2,y(0)=0}, y(x), numeric); nummego2 := proc(rkf45_x) ... end # A közelítő megoldás pl. x=1.3-ban kiolvasható az alábbi sorból: > nummego2(1.3); [x = 1.3, y(x) = .6610955532827162] # Tehát y(1.3) közelítőleg 0.6610955532827162 . # # A könnyebb használhatóság kedvéért csinálhatunk ebből egy "rendes" # közelítő függvényt: > ynum2 := x->rhs(nummego2(x)[2]); ynum2 := x -> rhs(nummego2(x)[2]) # Egy próba: > ynum2(1.3); .6610955532827162 # Ez ugyanaz, mint az előbb, ahogy vártuk. # # Ezzel az ynum2-vel nem végezhető el minden, amit megszoktunk, mert # csak # egyes helyeken értelmes, de pl. rajzoltatni tudjuk. Pl. > diff(ynum2(x), x); #hiba Error, (in nummego2) cannot evaluate boolean > plot('ynum2(x)', x=0..10); # Példa. # Rajzold ki az alábbi kezdeti érték feladat megoldásának grafikonját az # [1,20] intervallumon! # y'=(2-x)(x+y^2), y(1)=2 > diffegy:=diff(y(x), x)=(2-x)/(x+y(x)^2); d 2 - x diffegy := -- y(x) = --------- dx 2 x + y(x) > kezdfelt:=y(1)=2; kezdfelt := y(1) = 2 > dsolve( {diffegy,kezdfelt}, y(x)); # Mivel nem kaptunk függvényt megoldásként, ezért közelítö megoldásra # van szükség. # > nummegold:=dsolve( {diffegy,kezdfelt}, y(x), numeric); nummegold := proc(rkf45_x) ... end # Próba, hogy tényleg kaptunk-e valamit: > nummegold(1.5); [x = 1.5, y(x) = 2.069989241662746] > ynum:=x->rhs(nummegold(x)[2]); ynum := x -> rhs(nummegold(x)[2]) # Próba: > ynum(1.5); 2.069989241662746 > plot('ynum(x)', x=1..20); >