1 / 56

Riešenie rovníc v Matlabe

Riešenie rovníc v Matlabe. Lineárne numerické symbolické Nelineárne numerické symbolické Diferenciálne numerické symbolické. Kozák 2001-2004. Demonštračné príklady ku kurzu Matlabu 1999 - 2002 ver. 6. 5-7.1.

howard
Download Presentation

Riešenie rovníc v Matlabe

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Riešenie rovníc v Matlabe Lineárne numerické symbolické Nelineárnenumerické symbolické Diferenciálne numerické symbolické Kozák 2001-2004

  2. Demonštračné príklady ku kurzu Matlabu 1999-2002 ver.6.5-7.1 % Demo1 - pripravil Š.Kozak,KASR 2001-2004 %Numerické riešenie linearnych rovnic Ax=b, x=inv(A) A=[1 2 1;3 4 1;2 5 7] % vstup matice b=[1 2 3]' %vstup vektora ' znamena transponovanie x1=A\b % ries. systemu rovnic v x1 je riesenie echo on pause % prerus. behu, pokracovanie lub.klavesnicou format long % dlhy format (default je short) x2=inv(A)*b% iny sposob riesenia cez inverziu pause format short% kratky format (4 miesta za des.bodkou) E1=A*x1% skuska spravnosti pause E2=A*x2% skuska spravnosti

  3. % Symbolické a numerické riešenie lineárnych algebraických rovníc syms a11 a12 a13 a21 a22 a23 a31 a32 a33 b1 b2 b3; %deklaracia symbolickych premennych A1 = [a11 a12 ; a21 a22]; % zadanie matice vo vseobecnej forme (2x2) B1=[b1 b2];% zadanie vektora pravej strany XII = A1\B1'; % riesenie systemu rovnic x11 = XII(1) % vyber konkr.riesenia z celkoveho riesenia x12 = XII(2) pause A2=[a11 a12 a13;a21 a22 a23;a31 a32 a33] % zadanie matice vo vseobecnej forme (3x3) B2 = [b1 b2 b3]; a11=1;a12=2;a13=1;a21=3;a22=4;a23=1; a31=2;a32=5;a33=7;b1=1;b2=2;b3=3; XIII=A2\B2'; x21 = XIII(1) x22 = XIII(2) x33 = XIII(3) xs2=subs(XIII)% dosadenie (substitucia) do symbolickeho riesenia b=A*xs2 % skuska riesenia pause Ain=inv(A2) % symbolicka inverziA xs3=Ain*B2'% symbolicke riesenie xr=subs(xs3) % dosadenie (substitucia) do symbolickeho riesenia

  4. Symbolické riešenie rovníc Syntax: SOLVE Symbolické riešenie algebraických rovníc SOLVE('eqn1','eqn2',...,'eqnN') SOLVE('eqn1','eqn2',...,'eqnN','var1,var2,...,varN') SOLVE('eqn1','eqn2',...,'eqnN','var1','var2',...'varN')

  5. Symbolické riešenie sústavy lineárnych rovníc Priamy spôsob riešenia lin.rovníc (nezabudnúť príkaz syms-definovanie symbolických premenných) syms a11 a12 a21 a22 b1 b2 ; %deklaracia symbolickych premennych [x1,x2]=solve('a11*x1+a12*x2=b1','a21*x1+a21*x2=b2') %koniec symbolickeho riesenia, takto mozno riesit aj nelinearne systemy

  6. Symbolické riešenie nelineárnych rovníc % Kvadratická rovnica syms a b c x x = solve(a*x^2 + b*x + c) % sym.riesenie kvadr.rovnice x = solve('a*x^2 + b*x + c = 0') % Kubická rovnica syms a b c x x=solve('a*x^3+b*x^2+c*x+d')

  7. Symbolické riešenie transcedentálnych rovníc [x,y] = solve('sin(x+y)-exp(x)*y = 0','x^2-y = 2') [x1,x2]=solve('2*x1 - x2 - exp(-x1)',... '-x1 + 2*x2 - exp(-x2) ')

  8. Symbolické riešenie diferenciálnych rovníc y = dsolve('Dy = -a*y','y(0) = 1') % symb. ries. dif. rovnice 1.radu y = dsolve('D2y = -a^2*y', 'y(0) = 1, Dy(pi/a) = 0') % symb. ries. dif. rovnice 2. radu Počiatočné podmienky y = dsolve('D2y = sin(y)'); pretty(y)

  9. dsolve('eqn1','eqn2', ...) symb.ries. dif. rovnic Príklady: dsolve('Dx = -a*x') ans = exp(-a*t)*C1 x = dsolve('Dx = -a*x','x(0) = 1','s') x = exp(-a*s) y = dsolve('(Dy)^2 + y^2 = 1','y(0) = 0') y = [ sin(t)] [ -sin(t)] S = dsolve('Df = f + g','Dg = -f + g','f(0) = 1','g(0) = 2')

  10. S.f = exp(t)*cos(t)+2*exp(t)*sin(t) S.g = -exp(t)*sin(t)+2*exp(t)*cos(t) Y = dsolve('Dy = y^2*(1-y)') Warning: Explicit solution could not be found; implicit solution returned. Y = t+1/y-log(y)+log(-1+y)+C1=0 dsolve('Df = f + sin(t)', 'f(pi/2) = 0') dsolve('D2y = -a^2*y', 'y(0) = 1, Dy(pi/a) = 0') S = dsolve('Dx = y', 'Dy = -x', 'x(0)=0', 'y(0)=1') S = dsolve('Du=v, Dv=w, Dw=-u','u(0)=0, v(0)=0, w(0)=1') w = dsolve('D3w = -w','w(0)=1, Dw(0)=0, D2w(0)=0') y = dsolve('D2y = sin(y)'); pretty(y)

  11. Numerické riešenie nelinárnych rovníc FSOLVE riešenie nelineárnych rovníc metódou najmenších štvorcov F(X)=0kde F a X sú vektory alebo matice Syntax : X=FSOLVE(FUN,X0)% X0 vektor start. hodnot FUN system rovníc rovnica) X výsledný vektor alebo matica X=FSOLVE(FUN,X0,OPTIONS) alebo [X,FVAL,EXITFLAG]=FSOLVE(FUN,X0,...) If EXITFLAG je: > 0 riesenie konverguje FSOLVE converged to a solution X. = 0 potom max. počet iterácií bol dosiahnuty. < 0 riesenie nekonverguje

  12. 1. 2. 3. fun = inline('sin(3 *x)'); x = fsolve(fun,[1 4],optimset('Display','off')) x = 1.0472 4.1888 fun1 = inline('(x.^2+2.*x+1) ') x = fsolve(fun1,[0 4],optimset('Display','on')) Tretí spôsob riešenia NR: x = fsolve(@mojafun,[2 3 4],optimset('Display','iter')) kde MOJAFUNje MATLAB funkcia napr: function F = mojafun(x) F = sin(x);

  13. fsolve-príklad Pre dve NLR : 2x1-x2–e-x1 =0, -x1+2x2-e-x2=0 Poc.podmienky x0=[-5 -5]. Najlepsie je vytvorit si zvlast funkciu function F = mojafun(x) F = [2*x(1) - x(2) - exp(-x(1));       -x(1) + 2*x(2) - exp(-x(2))]; x0 = [-5; -5];       %Startovacie podmienky options=optimset('Display','iter'); % Zobrazovanie výsledkov (x-je riešenie, fval hodnota ucel.funkcie) Optimization terminated successfully: [x,fval] = fsolve(@mojafun,x0,options)

  14. Vypis riesenia : (x-je riešenie, fval hodnota ucel.funkcie) Optimization terminated successfully: [x,fval] = fsolve(@mojafun,x0,options) Directional Iteration Func-count Residual Step-size derivative 1 2 39.6268 1 -79.3 2 9 2.14377 1.44 11.5 3 15 0.00107228 1.06 -0.0155 4 21 6.6055e-012 1 -1.85e-008 Optimization terminated successfully: Search direction less than tolX x = 0.5671 0.5671 fval = 1.0e-008 * -0.2820 0.6111

  15. fzero Najdenie koreňa funkcie jednej premenneje Syntax x = fzero(fun,x0) x = fzero(fun,x0,options) x = fzero(fun,x0,options,P1,P2,...) [x,fval] = fzero(...) [x,fval,exitflag] = fzero(...) [x,fval,exitflag,output] = fzero(...)

  16. Príklady: X = FZERO(FUN,X0) FUN môže byť definovaná pomocou znaku @: X = fzero(@sin,3) Riešením je pi X = fzero(@sin, 3, optimset('disp','iter')) riešením je pi, (default tolerance a zobrazuje iterácie FUN môže byť zadávaná ako inline objekt: X = fzero(inline('sin(3*x)'),2); Limitations: X = fzero(inline('abs(x)+1'), 1) returns NaN since this function does not change sign any-where on the real axis (and does not have a zero as well).

  17. X = fzero(@tan,2) returns X near 1.5708 because the discontinuity of this function near the point X gives the appearance (numerically) that the function changes sign at X.

  18. Riešenie diferenciálnych rovníc [T,Y] = ODE45(ODEFUN,TSPAN,Y0) with TSPAN = [T0 TFINAL] [T,Y] = ODE45(ODEFUN,TSPAN,Y0,OPTIONS) ODE solvers:ODE23, ODE113, ODE15S, ODE23S, ODE23T, ODE23TB options handling: ODESET, ODEGET output functions: ODEPLOT, ODEPHAS2, ODEPHAS3, ODEPRINT evaluating solution: DEVAL ODE examples: RIGIDODE, BALLODE, ORBITODE Metódy integrácie :

  19. Príklady: [t,y]=ode45(@vanderpoldemo,[0 20],[2 0]) plot(t,y(:,1)) function dydt = vanderpoldemo(t,y,Mu) %VANDERPOLDEMO DefinujeVan der Polovu rovnicu. dydt = [y(2); Mu*(1-y(1)^2)*y(2)-y(1)]; [t,y] = ode23(@prstrana,[t0 tfinal],y0) function yp = prstrana(t,y) %vektor pravych stran. dif.rovnic %y1' = (1 - alpha*y2)*y1, alpha=0.01, beta=0.02 %y2' = (-1 + beta*y1)*y2 % yp=[y(1),y(2)]=… % zapisana ako yp‘=Ay, A=[1 - .01*y(2) 0;0 -1 + .02*y(1)] yp = diag([1 - .01*y(2), -1 + .02*y(1)])*y; % alebo iny zapis yp=A*y, ak dopredu definujeme A

  20. Demonštračné príklady integračné metódy% Model zásobníka vodyecho offglobal ro g A1 A2 Hh Qhro=1000 % [kg/m3]g=9.81 % [m/sek2]A1=pi/4; % [m2]A2=pi/400; % [m2]Hh=1; % [m]Qh=34.77 % [kg/sek]% Integracia systemu nelinearnych rovnic metódou ODE45% Vstupne parametre (nevyhnutne)% Pociatocny cas integracie t0t0=0;% Koncovy cas integracie tftf=300;H0=0.5 % Pociatocna podmienka[t,H]=ode45(@zasps2,[t0 tf],H0); % volanie procedury integraciefigure(1)plot(t,H);gridxlabel('Cas [sek]'),ylabel('H(m)') Q1(t) H Q2(t)

  21. title('Integracia systemu nelinearnych rovnic') %Vypocet vystupneho prietoku figure(2) k3=ro*sqrt(2*g)*A2; Q2=k3*sqrt(H); plot(t,Q2);grid xlabel('Cas [sek]'),ylabel('Q2 [kg/sek]') function Hder=zasps2(t,H)function Q=q(t) global ro g A1 A2 Hh Qh global Qh Q1=q1(t); Q=Qh; k1=-A2*sqrt(2*g)/A1; k2=1/ro/A1; Hder=k1*sqrt(H)+k2*Q1; function Hder=hder(t,H) global ro g A1 A2 Hh Qh Q=q(t) k1=-A2*sqrt(2*g)/A1 k2=1/ro/A1 Hder=k1*sqrt(H)+k2*Q function Q1=q1(t) global Qh k=1; w=0.05; %frekvencia sin.signalu DQ=k*sin(w*t); Q1=Qh+DQ;

  22. Príklady - riesenie a zapis diferencialnych rovnic % Uvazujme system dvoch difrenciálnych rovníc % y1' = (1 - alpha*y2)*y1 % y2' = (-1 + beta*y1)*y2 % Pociatocny cas integracie t0=0, koncový čas tf=15 % Default hodnota presnosti integracie je 1e-3 (0.1 percent). t0 = 0; tfinal = 15; y0 = [20 20]'; % Definicia poc.podmienok % Volanie procedury : % [t,y] = ode23(@prstrana,[t0 tfinal],y0); % funkcia prstrana je prava str.difer.rovnic, zapisuje sa % samostatne % [t,y]=ode45(…..) % Alebo ina metoda integracie ODE45, ODE113, ODE15S, % ODE23S, ODE23T, ODE23TB

  23. pause [t,y] = ode23(@prstrana,[t0 tfinal],y0); % vol.proc.integr. %vykreslenie trajektorii v casovej oblasti plot(t,y), title('Demo-priklad casova oblast') pause %vykreslenie trajektorii vo fazovej rovine plot(y(:,1),y(:,2)), title('Demo- priklad faz.rovina') pause function yp = prstrana(t,y) %vektor pravych stran. dif.rovnic %y1' = (1 - alpha*y2)*y1, alpha=0.01, beta=0.02 %y2' = (-1 + beta*y1)*y2 % yp=[y(1),y(2)]=… % zapisana ako yp‘=Ay, A=[1 - .01*y(2) 0;0 -1 + .02*y(1)] yp = diag([1 - .01*y(2), -1 + .02*y(1)])*y; % alebo iny zapis yp=A*y, ak dopredu definujeme A

  24. Iné priklady zápisu riesenia diferencialnych rovnic  Vhodny pre starsie aj nove verzie ODE % [T,Y] = ode45(@prstr1,[t0 tfinal],y0); %Diferencialna rovnica je druheho radu %napr.Van der Pol rovnica y'' + (y^2 - b)y' + y = 0 %v takomto pripade musime ju znizit na system rovnic prveho radu % y1' = y1(1 - y2^2) - y2 % y2' = y1 pause % Strike any key to continue. clc % To simulate the differential equation defined in VDPOL over the % interval 0 < t < 20, použitím ODE23:

  25. t0 = 0; tfinal = 15; y0 = [0 0.25]; % Define initial conditions. % [t,y] = ode23(@prstr1,t0,tfinal,y0); [t,y] = ode23(@prstr1,t0,tfinal,y0); % to je % najjednoduchsie volanie procedury integracie plot(t,y), title('Zobrazenie traj. v casovej oblasti'), pause plot(y(:,1),y(:,2)) title('Zobrazenie trajektorii vo fazovej rovine‘ pause clc % Ak chceme pouzit inu proced.integracie napr. ode 45 % [T,Y] = ode45(@vdpol,t0,tfinal,y0); %volanie proced.

  26. function yp = prstr1(t,y);  % prava strana je zapisana ako vektor jednotlivych dif.rovnic yp = [(y(1) .* (1 - y(2).^2) - y(2)); y(1)]; Example 1. Trojrozmerný dynamický systém Zápis pravých stran: function dy = rigid(t,y) dy = zeros(3,1); % a column vector dy(1) = y(2) * y(3); dy(2) = -y(1) * y(3); dy(3) = -0.51 * y(1) * y(2); %Nastavenie presnosti riešenia options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-5]); [t,y] = ode45(@rigid,[0 12],[0 1 1],options); Plotting the columns of the returned array Y versus T shows the solution: plot(t,y(:,1),'-',t,y(:,2),'-.',t,y(:,3),'.')

  27. Rozšírenia neprednášať

  28. [x,y]=solve('x^2*y^2 - 2*x - 1 = 0','x^2 - y^2 - 1 = 0') % riesenie systemu 2NR

  29. [x1,x2]=solve('a11*x1+a12*x2=b1','a21*x1+a21*x2=b2') % priamy sposob riesenia lin.rovnic %koniec symbolickeho riesenia, takto mozno riesit aj nelinearne systemy syms a b c x x = solve(a*x^2 + b*x + c)% symbol. riesenie kvadr.rovnice x = solve('a*x^2 + b*x + c = 0') [x,y]=solve('x^2*y^2 - 2*x - 1 = 0','x^2 - y^2 - 1 = 0') % riesenie % systemu [x,y] = solve('sin(x+y)-exp(x)*y = 0','x^2-y = 2') y = dsolve('Dy = -a*y','y(0) = 1') % symb. ries. dif. rovnic y = dsolve('D2y = -a^2*y', 'y(0) = 1, Dy(pi/a) = 0') y = dsolve('D2y = sin(y)'); pretty(y)

  30. Príklad : Riešenie systému nelinearnych rovnic funkcia echo on % Demo2 - pripravil Š.Kozak,KASR 1999 %Klasické riešenie nelinearnych rovnic x=fsolve('x^2+2*x+1', 0) pause x=fsolve(@(x^2+2*x+1), 0) % symbolicke riesenie pause solve('a*x^3+b*x^2+c*x+d') solve('a*x^3+b*x^2+c*x+d') solve('x^3+3*x^2+3*x+1')

  31. fsolve Pre dve NLS : 2x1-x2–e-x1 =0, -x1+2x2-e-x2=0 Poc.podmienky x0=[-5 -5]. Najlepsie je vytvorit si zvlast funkciu function F = mojafun(x) F = [2*x(1) - x(2) - exp(-x(1));       -x(1) + 2*x(2) - exp(-x(2))];  x0 = [-5; -5];       %Startovacie podmienky options=optimset('Display','iter'); % Zobrazovanie výsledkov (x-je riešenie, fval hodnota ucel.funkcie) Optimization terminated successfully: [x,fval] = fsolve('mojafun',x0,options)%ver.5.3

  32. Vypis riesenia : (x-je riešenie, fval hodnota ucel.funkcie) Optimization terminated successfully: [x,fval] = fsolve('mojafun',x0,options) Directional Iteration Func-count Residual Step-size derivative 1 2 39.6268 1 -79.3 2 9 2.14377 1.44 11.5 3 15 0.00107228 1.06 -0.0155 4 21 6.6055e-012 1 -1.85e-008 Optimization terminated successfully: Search direction less than tolX x = 0.5671 0.5671  fval = 1.0e-008 * -0.2820 0.6111

  33. V starších verziach použijeme priamo zápis riešenia F = '[2*x(1) - x(2) - exp(-x(1));-x(1)+2*x(2)-exp(-x(2))]' [x,fval]=fsolve(F,[-0.5 -0.5]) Riešenie je možné získať priamo : [x,fval]=fsolve('[2*x(1) - x(2) - exp(-x(1));-x(1)+2*x(2)-exp(-x(2))]',[-0.5 -0.5]) x = 0.5671 0.5671 fval = 1.0e-009 * -0.2460 -0.2460

  34. X=FSOLVE(FUN,X0) • Priklady .6.1 • FUN moze byt zavedena cez @: • x = fsolve(@myfun,[2 3 4],optimset('Display','iter')) • kde MYFUN je MATLABovska funkcia: • function F = myfun(x) • F = sin(x); • FUN moze byt zavedena cez inline prikaz: • fun = inline('sin(3*x)'); • x = fsolve(fun,[1 4],optimset('Display','off'))

  35. Finl=inline('[2*x(1) - x(2) - exp(-x(1))]; -x(1) + 2*x(2) - exp(-x(2))]') Finl = Inline function: Finl(x) = [2*x(1) - x(2) - exp(-x(1))]; -x(1) + 2*x(2) - exp(-x(2))] x0=[-5 -5]      %Startovacie podmienky options=optimset('Display','iter'); % Zobrazovanie výsledkov  [x,fval] = fsolve(Finl,x0,options) % volanie opt.funkcie Finl=inline('2*x(1) - x(2) - exp(-x(1), -x(1) + 2*x(2) - exp(-x(2)') Finl = Inline function: Finl(x) = [2*x(1) - x(2) - exp(-x(1))]; -x(1) + 2*x(2) - exp(-x(2))]

  36. Finl=inline('2*x(1) - x(2) - exp(-x(1)); -x(1) + 2*x(2) - exp(-x(2))') • Finl = • Inline function: • Finl(x) = [2*x(1) - x(2) - exp(-x(1))]; -x(1) + 2*x(2) - exp(-x(2))]

  37. fzero Find zero of a function of one variable Syntax x = fzero(fun,x0) x = fzero(fun,x0,options) x = fzero(fun,x0,options,P1,P2,...) [x,fval] = fzero(...) [x,fval,exitflag] = fzero(...) [x,fval,exitflag,output] = fzero(...)

  38. X = FZERO(FUN,X0)Examples FUN can be specified using @: X = fzero(@sin,3) returns pi. X = fzero(@sin, 3, optimset('disp','iter')) returns pi, uses the default tolerance and displays iteration information. FUN can also be an inline object: X = fzero(inline('sin(3*x)'),2); Limitations X = fzero(inline('abs(x)+1'), 1) returns NaN since this function does not change sign anywhere on the real axis (and does not have a zero as well).

  39. X = fzero(@tan,2) returns X near 1.5708 because the discontinuity of this function near the point X gives the appearance (numerically) that the function changes sign at X.

  40. % Find a zero of humps(x) near 1 with FZERO z = fzero(@humps,1,optimset('Display','off')) fplot(@humps,[0,2]) hold on; plot(z,0,'r*'); hold off

  41. Q1(t) Demonštračné príklady integračné metódy% Model zásobníka vodyecho offglobal ro g A1 A2 Hh Qhro=1000 % [kg/m3]g=9.81 % [m/sek2]A1=pi/4; % [m2]A2=pi/400; % [m2]Hh=1; % [m]Qh=34.77 % [kg/sek]% Integracia systemu nelinearnych rovnic metódou ODE45% Vstupne parametre (nevyhnutne)% Pociatocny cas integracie t0t0=0;% Koncovy cas integracie tftf=300;H0=0.5 % Pociatocna podmienka[t,H]=ode45('zasps2',[t0 tf],H0); % volanie procedury integraciefigure(1)plot(t,H);gridxlabel('Cas [sek]'),ylabel('H(m)') H Q2(t)

  42. title('Integracia systemu nelinearnych rovnic') %Vypocet vystupneho prietoku figure(2) k3=ro*sqrt(2*g)*A2; Q2=k3*sqrt(H); plot(t,Q2);grid xlabel('Cas [sek]'),ylabel('Q2 [kg/sek]') function Hder=zasps2(t,H)function Q=q(t) global ro g A1 A2 Hh Qh global Qh Q1=q1(t); Q=Qh; k1=-A2*sqrt(2*g)/A1; k2=1/ro/A1; Hder=k1*sqrt(H)+k2*Q1; function Hder=hder(t,H) global ro g A1 A2 Hh Qh Q=q(t) k1=-A2*sqrt(2*g)/A1 k2=1/ro/A1 Hder=k1*sqrt(H)+k2*Q function Q1=q1(t) global Qh k=1; w=0.05; %frekvencia sin.signalu DQ=k*sin(w*t); Q1=Qh+DQ;

  43. Riesenie a zapis diferencialnych rovnic % Uvazujme system dvoch difrenciálnych rovníc % y1' = (1 - alpha*y2)*y1 % y2' = (-1 + beta*y1)*y2 % % Pociatocny cas integracie t0=0, koncový čas tf=15 % Default hodnota presnosti integracie je 1e-3 (0.1 percent). t0 = 0; tfinal = 15; y0 = [20 20]'; % Definicia poc.podmienok % Volanie procedury : % [t,y] = ode23(@prstrana,[t0 tfinal],y0); % funkcia prstrana je prava str.difer.rovnic, zapisuje sa % samostatne % [t,y]=ode45(…..) % Alebo ina metoda integracie ODE45, ODE113, ODE15S, % ODE23S, ODE23T, ODE23TB

  44. pause [t,y] = ode23(@prstrana,[t0 tfinal],y0); % takto volame % proc.integ. %vykreslenie trajektorii v casovej oblasti plot(t,y), title('Demo-priklad casova oblast') pause %vykreslenie trajektorii vo fazovej rovine plot(y(:,1),y(:,2)), title('Demo- priklad fazova rovina'),pause function yp = prstrana(t,y) %vektor pravych stran. dif.rovnic %y1' = (1 - alpha*y2)*y1, alpha=0.01, beta=0.02 %y2' = (-1 + beta*y1)*y2 % yp=[y(1),y(2)]=… % zapisana ako yp‘=Ay, A=[1 - .01*y(2) 0;0 -1 + .02*y(1)] yp = diag([1 - .01*y(2), -1 + .02*y(1)])*y; % alebo iny zapis yp=A*y, ak dopredu definujeme A

  45. Iné priklady zápisu riesenia diferencialnych rovnic  Vhodny pre starsie aj nove verzie ODE % [T,Y] = ode45('prstr1',[t0 tfinal],y0); %Diferencialna rovnica je druheho radu %napr.Van der Pol rovnica y'' + (y^2 - b)y' + y = 0 %v takomto pripade musime ju znizit na system rovnic prveho radu % y1' = y1(1 - y2^2) - y2 % y2' = y1 pause % Strike any key to continue. clc % To simulate the differential equation defined in VDPOL over the % interval 0 < t < 20,použitímODE23:

  46. t0 = 0; tfinal = 15; y0 = [0 0.25]; % Define initial conditions. % [t,y] = ode23('prstr1',t0,tfinal,y0); [t,y] = ode23('prstr1',t0,tfinal,y0); % to je % najjednoduchsie volanie procedury integracie plot(t,y), title('Zobrazenie traj. v casovej oblasti'), pause plot(y(:,1),y(:,2)) title('Zobrazenie trajektorii vo fazovej rovine‘ pause clc %Ak chceme pouzit inu proceduru integracie napr. ode 45 % [T,Y] = ode45('vdpol',t0,tfinal,y0); _ volanie procedury

  47. function yp = prstr1(t,y);  % prava strana je zapisana ako vektor jednotlivych dif.rovnic yp = [(y(1) .* (1 - y(2).^2) - y(2)); y(1)]; Example 1. Trojrozmerný dynamický systém Zápis pravých stran: function dy = rigid(t,y) dy = zeros(3,1); % a column vector dy(1) = y(2) * y(3); dy(2) = -y(1) * y(3); dy(3) = -0.51 * y(1) * y(2); %Nastavenie presnosti riešenia options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-5]); [t,y] = ode45('rigid',[0 12],[0 1 1],options); Plotting the columns of the returned array Y versus T shows the solution: plot(T,Y(:,1),'-',T,Y(:,2),'-.',T,Y(:,3),'.')

  48. Ako je to v novej verzii Matlabu ? Example 2 verziaMatlab 6.0 [t,y]=ode45(@vdp1,[0 20],[2 0]); plot(t,y(:,1)); solves the system y' = vdp1(t,y), using the default relative error tolerance 1e-3 and the default absolute tolerance of 1e-6 for each component, and plots the first component of the solution. ODE solvers: ODE23, ODE113, ODE15S, ODE23S, ODE23T, ODE23TB options handling: ODESET, ODEGET output functions: ODEPLOT, ODEPHAS2, ODEPHAS3, ODEPRINT ODE examples: RIGIDODE, BALLODE, ORBITODE

More Related