ˇ
Rešení
diferenciálních rovnic v MATLABu
Základy algoritmizace a programování
Pˇrednáška
23. listopadu 2011
ZAPG
Co ˇrešíme
Obyˇcejné diferenciální rovnice
prvního ˇrádu:
separovatelné
lineární
exaktní
druhého ˇrádu, lineární:
s konstantními koeficienty
ˇ
s promennými
koeficienty
soustavy
lineární: X = AX + B
nelineární
ZAPG
Obyˇcejné diferenciální rovnice
ˇ
Rešení
Cauchyovy úlohy
y 0 = f (t, y (t)), y (t0 ) = t0 ,
kde t je skalár, y (t) neznámý vektor hodnot, f (t, y ) zadaná
funkce (resp. vektor hodnot).
ZAPG
ˇ
ˇ
Rešení
v symbolických promenných
Obyˇcejná diferenciální rovnice 1. ˇrádu, napˇr.
y0 +
1
y
=
x
x +3
y (−2) = 4
1
rovnici upravíme na tvar : y 0 = f (x, y ).
2
help dsolve
3
funkce dsolve(’Dy =f(x)’,’y(x0)=y0’, ’x’)
y=dsolve(’Dy=1/(x+3)-y/x’,’y(-2)=4’,’x’)
y je výraz, který odpovídá ˇrešení rovnice a do kterého lze
dosadit, napˇr. vektor hodnot
4
zobrazení pˇresného ˇrešení
xx=-2:0.1:-0.1;
yy=subs(y,xx);
plot(xx,yy)
nebo
ezplot(y)
ZAPG
v matlabu
y =
(x - 3*log(x + 3))/x - 6/x
latex(y)
ans=
x − 3 ln(x + 3) 6
−
x
x
ZAPG
Lineární diferenciální rovnice 2. ˇrádu s konstantními koeficienty,
y 00 − 9y = 5e2t ,
y (0) = 0, y 0 (0) = 3
y = dsolve(’D2y = 5*exp(2*t)+9*y’,’y(0)=0, Dy(0)=3’,’t’)
ezplot(y)
ZAPG
v matlabu
y =
(4*exp(3*t))/3 - 1/(3*exp(3*t)) - exp(2*t)
» latex(y)
ans =
1
4 e3 t
− 3 t − e2 t
3
3e
ZAPG
ˇ
Lineární diferenciální rovnice 2. ˇrádu s promennými
koeficienty,
y 00 + xy =
1
,
2−x
y (0) = 0.25, y 0 (0) = 0
y = dsolve(’D2y = 1/(2-x)-x*y’, ’y(0)=0.25, Dy(0)=0’, ’x’)
Ne vždy je ˇrešení nalezeno . . .
ZAPG
Lineární soustavy,
X˙ =
1
1
4 −2
X,
X
0
5
dsolve(’Dy = x+y, Dy = 4*x -2y’, ’x(0)=0, y(0)=5’, ’t’)
ZAPG
Numerické ˇrešení
ˇ
Napíšeme funkci, která poˇcítá hodnoty f a použijeme nekterou
z matlabovských funkcí, které potˇrebují:
jméno funkce,
rozsah hodnot t (od t0 do tN ) a
poˇcáteˇcní hodnotu y0 .
ˇ matlabovská funkce: ode45.
Nejpoužívanejší
(help ode45)
ZAPG
Použití ode45
y 0 = −y (t) − 5e−t sin 5t, y (0) = 1, pro 0 ≤ t ≤ 3.
Vytvoˇríme funkci prava (v souboru prava.m)
function dy = prava(t,y)
dy = -y-5*exp(-t)*sin(5*t);
a použijeme ji:
interval_t = [0 3]; y0 = 1;
[t, y ] = ode45(@prava, interval_t, y0);
ˇ
v promenné
t jsou body z intervalu < 0, 3 >, ve kterých jsou
urˇceny hodnoty pˇribližného rˇešení y(i).Vnitˇrní body vybírá
ˇ ˇrešení.
funkce ode45, tím menší vzdálenost, cˇ ím více se mení
Výsledek mužeme
˚
graficky zobrazit:
plot(t,y,’*- -’);
xlabel t, ylabel y(t)
Pro uvedený pˇríklad známe pˇresné ˇrešení: y (t) = e−x cos 5t,
proto mužeme
˚
urˇcit chybu numerického ˇrešení jako
max(abs(y-exp(-t).* cos(5*t))) = 2.8991 e -04 a porovnat grafy.
ZAPG
Pokud je zadáno více než 2 hodnoty (pro t) výpoˇcet je
ˇ
ˇ
proveden pouze v techto
hodnotách, a funkce delení
intervalu
neprovádí (respektuje zadané).
Napˇr.
tspan2 = 0:4
budou y hodnoty vypoˇcteny pro t=0,1,2,3,4
Lze zadat i obrácené uspoˇrádání hodnot (záporný krok), napˇr.
tspan3 = [0 -0.5 -1]
ZAPG
ODR vyššího ˇrádu
ˇ
Rešíme
pˇrevedením rovnice vyššího ˇrádu na soustavu
diferenciálních rovnic prvního rˇádu.
Napˇr. rovnice kyvadla:
y 00 = − sin y
Pˇrevedeme na soustavu 2 rovnic: y1 (t) = y (t) a y2 (t) = y 0 (t):
y10 (t) = y2 (t),
y20 (t) = − sin y1 (t)
Pro použití ode45 vytvoˇríme funkci pravé strany:
function dy = kyvadlo(t,y)
dy =[y(2); -sin(y(1))];
ZAPG
Kyvadlo
Výpoˇcet provedeme pro 0 ≤ t ≤ 10 s ruznými
˚
poˇcáteˇcními
podmínkami.
Návratovými hodnotami ode45 bude matice, která má v
každém ˇrádku t(i), y1(t), y2(t).
tspan = [0 10 ];
yA0 = [1; 1]; yB0 = [-5; 2 ]; yC0=[5; -2];
[ta ya] = ode45(@kyvadlo, tspan, yA0);
[tb yb] = ode45(@kyvadlo, tspan, yB0);
[tc yc] = ode45(@kyvadlo, tspan, yC0);
ZAPG
Zobrazení ˇrešení
Zobrazení fázových trajektorií v rovineˇ – v osách y1(t), y2(t) –
použijeme vygenerované sloupce y(:,1), y(:,2).
ˇ
Vektorové pole smerových
vektoru˚ [y2, -sin y1 ] zobrazí funkce
quiver.
[y1,y2] = meshgrid (-5:0.5:5, -3:0.5:3);
Dy1Dt = y2; Dy2Dt = -sin(y1);
quiver(y1,y2,Dy1Dt,Dy2Dt)
hold on plot(ya(:,1),ya(:,2))
plot(yb(:,1),yb(:,2))
plot(yc(:,1),yc(:,2))
axis equal, axis([-5 5 -3 3])
xlabel y_1(t), ylabel y_2(t), hold off
ZAPG
Kyvadlo
ZAPG
Soustavy autonomních rovnic
1
2
3
4
bod rovnováhy : ohnisko
1 −1
˙
X =
X,
1
1
bod rovnováhy : sedlo
1 −1
˙
X =
X,
−4 −2
bod rovnováhy : uzel
−1
0
˙
X =
X,
3 −2
bod rovnováhy : stˇred
0 1
˙
X =
X,
−4 0
ZAPG
X (0) =
0
5
1
0
X (0) =
X (0) =
X (0) =
1
2
1
0
Autonomní soustava: bod rovnováhy OHNISKO
tspan=[0,-5];
ybzero=[1;2];
[tb,yb]=ode45(@pr31,tspan,ybzero);
[y1,y2] = meshgrid(-1:0.2:2,-1:0.2:2);
Dy1Dt = y1-y2; Dy2Dt =y1+y2;
quiver(y1,y2,Dy1Dt,Dy2Dt);
hold on
plot(yb(:,1),yb(:,2))
axis equal
axis([-1,2,-1,2])
xlabel x(t), ylabel y(t), hold off
function yprime = pr31 (t,y)
yprime =[y(1)-y(2);y(1)+y(2)];
ZAPG
Ohnisko
ZAPG
Autonomní soustava: bod rovnováhy SEDLO
tspan=[0,1];
ybzero=[0;1]; yczero=[0; -1]; yazero=[1;0]; ydzero=[-1;0];
[tb,yb]=ode45(@pr19,tspan,ybzero);
[tc,yc]=ode45(@pr19,tspan,yczero);
[tb,ya]=ode45(@pr19,[0,0.3],yazero);
[tb,yd]=ode45(@pr19,[0,0.3],ydzero);
[y1,y2] = meshgrid(-1.5:0.2:1.5,-1.5:0.2:1.5);
Dy1Dt = y1-y2; Dy2Dt =-4*y1-2*y2;
quiver(y1,y2,Dy1Dt,Dy2Dt);
hold on
plot(yb(:,1),yb(:,2), yc(:,1),yc(:,2),ya(:,1),ya(:,2),yd(:,1),yd(:,2));
axis equal
xlabel y_1(t), ylabel y_2(t), hold off
function yprime = pr19 (t,y)
yprime =[y(1)-y(2);-4*y(1)-2*y(2)];
ZAPG
Sedlo
ZAPG
Autonomní soustava: bod rovnováhy UZEL
tspan =[0, 20];
yazero = [1; 0];
[ta,ya] = ode45(@pr12, tspan, yazero);
[y1,y2] = meshgrid(-1:0.2:1,-1:0.2:1);
Dy1Dt = -y2; Dy2Dt =3*y1-2*y2;
quiver(y1,y2,Dy1Dt,Dy2Dt);
hold on
plot(ya(:,1),ya(:,2),’r’)
axis equal
xlabel y_1(t), ylabel y_2(t), hold off
function yprime = pr12 (t,y)
yprime =[-y(2);3*y(1)-2*y(2)];
ZAPG
Uzel
ZAPG
ˇ
Autonomní soustava: bod rovnováhy STRED
tspan = [0,pi];
yazero = [1;0]; ybzero = [0; 3]; yczero = [1; 1];
[ta,ya] = ode45(@pr28, tspan, yazero);
[tb,yb] = ode45(@, tspan, ybzero);
[tc,yc] = ode45(@, tspan, yczero);
[y1,y2] = meshgrid(-1.5:0.3:1.5,-3:0.3:3);
Dy1Dt = y2; Dy2Dt = -4*y1;
quiver(y1,y2,Dy1Dt,Dy2Dt);
hold on
plot(ya(:,1),ya(:,2),yb(:,1),yb(:,2),yc(:,1),yc(:,2))
axis equal, axis([-2,2,-4,4])
xlabel y_1(t), ylabel y_2(t), hold off
function yprime = pr28 (t,y)
yprime =[y(2);-4*y(1)];
ZAPG
Stˇred
ZAPG
Numerická integrace v MATLABu
Numerická integrace – integrace v kvadraturách – pˇribližný
Rb
výpoˇcet a f (x)dx
funkce quad – realizuje Simpsonovu metodu
ˇ metodu
funkce quad1 – realizuje pˇresnejší
(Gauss–Lobato, Kronrod)
ˇ
funkce trapz – realizuje lichobežníkovou
metodu
funkce dblquad výpoˇcet dvojného integrálu
ZAPG
Použití funkcí
trapz :
parametry – vektory x– a y – souˇradnic, napˇr.
R 2π sin2 (x)
√
dx
0
2
1+cos (x)
>>x = linspace(0, 2*pi, 10); ... nebo x = 0: 2*pi/10: 10;
>>y = sin(x).ˆ2./sqrt(1+cos(x).ˆ2);
>>trapz(x,y)
ans =
2.8478
ZAPG
Použití
quad, quad1, dblquad
parametry: funkce, a, b, pˇresnost
funkce – musí mít parametr – vektor a vracet vektor funkˇcních
hodnot
(pro dblquad 2 parametry: vektor x, skalár y, vrací vektor
q = quad (funkce, a, b, presnost)
q = quad1(inline(’cos(x.ˆ2)’), t(i), t(i+1), 1e-3);
ZAPG
Pˇríklady
R4
x ln xdx
function f = xlnx(x)
f=x.*log(x);
2
quad(@xlnx,2,4)
>> ans =
6.7041
ZAPG
Pˇríklady
Rt
Rt
x(t) = 0 cos(u 2 )du, y (t) = 0 sin(u 2 )du
n=1000;
x=zeros(1,n); y=x;
t=linspace(0,4*pi,n+1);
for i=1:n
x(i)=quad(inline(’cos(x.ˆ2)’), t(i),t(i+1),1e-3);
y(i)=quad(inline(’sin(x.ˆ2)’), t(i),t(i+1),1e-3);
end
ZAPG
Obrázek
ZAPG
pˇríklady
R6R1
4
0
(y 2 ex + x cos(x))dxdy
function vysl = fxy(x,y)
vysl
=
yˆ2*exp(x)+x*cos(y);
>>dblquad(@fxy,0,1,4,6)
ans =
87.2983
ZAPG
Download

Rešení diferenciálních rovnic v MATLABu