Obyčejné diferenciální rovnice 1. řádu
Počáteční (Cauchyova) úloha pro jednu rovnici: y ' = f(x,y)
Znázornění všech řešení rovnice
Příklad: chceme znázornit všechna řešení rovnice y ' = y / x2 na obdélníkové oblasti [ 0.8 , 2.5 ] x [ 0.8 , 3.8 ] :
>> h = 0.1; % krok site >> xh = 0.8 : h : 2.5 ; % xh je interval [ 0.8 , 2.5 ] rozdeleny s krokem h >> yh = 0.8 : h : 3.8 ; % yh je interval [ 0.8 , 3.8 ] rozdeleny s krokem h >> [xx,yy] = meshgrid(xh,yh); % vytvoreni site na obdelniku [ 0.8 , 2.5 ] x [ 0.8 , 3.8 ] >> px = ones(length(yh),length(xh)); % px = 1 v kazdem bode [x,y] site >> py = yy./xx.^2; % py = f(x,y) v kazdem bode [x,y] site >> pp = sqrt(px.*px + py.*py); % (ppx, ppy) bude normovany vektor (px, py) >> ppx = px./pp; >> ppy = py./pp; >> quiver(xh, yh, px, py) % v kazdem bode [x,y] se nakresli vektor (ppx, ppy) jako sipka >> hold on % pokracovani ve stejnem obrazku >> plot(xh, 2*exp(1-1./xh),'LineWidth',2,'Color','r') % presne reseni >> axis([0.5, 2.6, 1.5, 3.6]) % nastaveni velikosti
Numerické řešení Cauchyovy úlohy y ' = f(x,y), y(x0) = y0
Příklad: přidáme do předchozího obrázku řešení rovnice y ' = y / x2 s počáteční podmínkou y(1) = 2. Rovnici řešíme numericky s krokem h = 0.2 buď Eulerovou nebo Collatzovou metodou.
>> h = 0.2; % krok >> x = 1; y = 2; % pocatecni podminka - vychozi bod >> hold on % chceme kreslit do stejneho obrazku >> plot(x, y, 'r*') % nakreslime vychozi bod jako cervenou hvezdicku
Eulerova metoda
Opakujeme následující 4 kroky:
>> f = y / x^2; % f ... smernice tecny v bode [x, y] >> x = x + h; % nova hodnota x >> y = y + f * h; % nova hodnota y (na primce se smernici f) >> plot(x, y, 'r*') % nakresleni noveho bodu
Collatzova metoda
Opakujeme následující kroky:
>> f = y / x^2; % f ... smernice tecny v bode [x, y] >> xp = x + 0.5 * h; % pomocny bod [xp, yp] >> yp = y + 0.5 * f * h; % (urceny Eulerovou metodou s polovicnim krokem) >> fp = yp / xp^2; % fp ... smernice tecny v bode [xp, yp] >> x = x + h; % nova hodnota x >> y = y + fp * h; % nova hodnota y (na primce se smernici fp) >> plot(x, y, 'r*') % nakresleni noveho bodu
Počáteční (Cauchyova) úloha pro soustavu 2 rovnic
Numerické řešení Cauchyovy úlohy Y ' = F(x,Y), Y(x0) = Y0
Příklad: numerické řešení rovnice Y ' = [ 2x - y2 + ln(x-1) + 1 , x (4-y1)1/2 - 1] T s počáteční podmínkou Y(2) = [ 0 , -1 ] T s krokem h = 0.1.
>> h = 0.1; % krok >> x = 2; Y = [ 0 ; -1 ]; % pocatecni podminky
Eulerova metoda
Opakujeme následující 3 kroky:
>> F = [ 2*x-Y(2)+log(x-1)+1 ; x*sqrt(4-Y(1))-1 ]; >> x = x + h; % nova hodnota x >> Y = Y + h * F; % nova hodnota Y
Collatzova metoda
Opakujeme následující kroky:
>> F = [ 2*x-Y(2)+log(x-1)+1 ; x*sqrt(4-Y(1))-1 ]; >> xp = x + 0.5 * h; % pomocny bod [xp, Yp] >> Yp = Y + 0.5 * h * F; >> Fp = [ 2*xp-Yp(2)+log(xp-1)+1 ; xp*sqrt(4-Yp(1))-1 ]; >> x = x + h; % nova hodnota x >> Y = Y + h * Fp; % nova hodnota Y
Převedení rovnice 2. řádu na soustavu 2 rovnic prvního řádu
Příklad - tlumený oscilátor: y''(x) + 2 y'(x) + y(x) = e-x, y(0) = 2 , y'(0) = -4 .
Označíme y1 = y, y2 = y' (tj. máme teď dvě nezávislé funkce: y1 představuje výchylku a y2 rychlost).
Platí y1' = y2, y2' = e-x - y1 - 2 y2, máme tedy soustavu diferenciálních rovnic
Y ' = [ y2 , e-x - y1 - 2 y2] T s počáteční podmínkou Y(0) = [ 2 , -4 ] T
Soustavu můžeme řešit jako výše, např. Eulerovou metodou s krokem h = 0.1 (v praxi spíš 0.01):
>> h = 0.1; % krok >> x = 0; Y = [ 2 ; -4 ]; % pocatecni podminky
Opakujeme následující 3 kroky:
>> F = [ Y(2) ; exp(-x)-Y(1)-2*Y(2) ]; >> x = x + h; % nova hodnota x >> Y = Y + h * F; % nova hodnota Y
Autonomní systém
Příklad - rovnice "dravec-kořist":
Y ' = [ a y1 - b y1 y2 ,
- c y2 + d y1 y2 ]
znázornění stavové (fázové) roviny na čtverci [ 0 , 40 ] x [ 0 , 40 ] :
>> a = 1.5; b = 0.2; c= 1; d = 0.1; % volba parametru >> h = 2; % krok site >> y1h = 0 : h : 30; % y1h je interval [ 0 , 40 ] rozdeleny s krokem h >> y2h = 0 : h : 40; % y2h je interval [ 0 , 40 ] rozdeleny s krokem h >> [y1,y2] = meshgrid(y1h,y2h); % vytvoreni site na obdelniku [ 0 , 40 ] x [ 0 , 40 ] >> py1 = a*y1 - b*y1.*y2; >> py2 = -c*y2 + d*y1.*y2; >> quiver(y1h, y2h, py1, py2, 2)
Eulerova metoda
>> h=0.1; >> Y = [15; 15]; % vychozi hodnota Y >> hold on >> plot(Y(1), Y(2), 'r*')
Opakujeme následující 3 kroky:
>> F = [ a*Y(1) - b*Y(1)*Y(2) ; -c*Y(2) + d*Y(1)*Y(2) ]; >> Y = Y + h * F; % nova hodnota Y >> plot(Y(1), Y(2), 'r*')
Collatzova metoda
>> h=0.1; >> Y = [15; 15]; % vychozi hodnota Y >> hold on >> plot(Y(1), Y(2), 'r*')
Opakujeme následující kroky:
>> F = [ a*Y(1) - b*Y(1)*Y(2) ; -c*Y(2) + d*Y(1)*Y(2) ]; >> Yp = Y + 0.5* h * F; % pomocny bod Yp >> Fp = [ a*Yp(1) - b*Yp(1)*Yp(2) ; -c*Yp(2) + d*Yp(1)*Yp(2) ]; >> Y = Y + h * Fp; % nova hodnota Y >> plot(Y(1), Y(2), 'r*')