dx dt = f x,t ( ) t

Size: px

Start display at page:

Download "dx dt = f x,t ( ) t"

ああすいまいだ
5 years ago
Views:

1 MATLAB Runge-Kutta

2 dx dt = f x,t ( ) t

3 dx( t) dt = lim Δt x( t + Δt) x( t) Δt t = nδt n = 0,1,2,3,4,5, x = x( nδt) n Δt dx ( ) dt = f x,t x n +1 x n Δt = f ( x,nδt) n

4 1 x n = x 0 n = 0 2 x = x + Δtf ( x,nδt) n +1 n n n = n +1 t = nδt Δt

5 Runge-Kutta x = x + Δtf x,t + Δt n +1 n 2 t = nδt x = x n + Δt ( ) 2 f x n,t

6 Euler.m Runge-Kutta RK2.m

7 Euler.m x0 x

8 Runge-Kutta RK2.m ( ) x0 x Runge-Kutta

9 MATLAB Runge-Kutta

10 τ dx ( t ) dt = x( t) + u( t) x L{ u( t) } = U( s) L{ x( t) } = X( s) ( ) = X ( s ) U( s) = 1 G s 1+ τs

11 U( s) = 1 s X( s) = 1 1+ τs 1 s = 1 s 1 s +1 τ x( t) =1 exp( t τ)

13 MATLAB Runge-Kutta

14 Euler.m Runge-Kutta RK2.m FOL.m FOL_main.m

15 FOL.m dx( t) dt = x ( t ) + u( t) τ

16 FOL_main.m function main tau = 2; % In = 1; % dt = 0.3; % T = 10; % Runge- Kutta xinit =[0]; % para = [tau,in]; % %@FOL [t1,x] = Euler(@FOL,[0,T],xinit,dt,para); % xe = x(1,:); xmin = min(xe)-0.1; xmax = max(xe)+0.1; %@FOL [t2,x] = RK2(@FOL,[0,T],xinit,dt,para); %2 Runge-Kutta xr2 = x(1,:); t3 = 0:dt:T; % xas = In.*(1-exp(-t3./tau)); % plot(t3,xas, -,t1,xe, -.,t2,xr2, : ); % axis([0 T xmin xmax]); end

17 MATLAB Runge-Kutta

18 Runge-Kutta RK2.m FOL_FR.m FOL_FR_main.m

19 FOL_FR.m function x = FOL_FR(t,y,para) x = (-y(1) + para(2).*sin(para(3).*t))./para(1); end dx( t) dt = x ( t ) + Asin( ωt) τ

20 FOL_FR_main. m function main tau = 1; % In = 1; % omega = 10; dt = 0.02; % T = 20; % xinit =[0]; % para = [tau,in,omega]; % %@FOL_FR [t,x] = RK2(@FOL_FR,[0,T],xinit,dt,para); %2 Runge-Kutta xr = x(1,:); xmin = min(xr)-0.1; xmax = max(xr)+0.1; plot(t,xr); % axis([0 T xmin xmax]);

21 dx( t) dt = x ( t ) + sin( ωt) τ ( ) x( t) = G( jω) sin ωt + tan 1 G( jω) G( jω) = 1 1+ jτω

23 TF_FOL.m FOL_boad_main.m

24 TF_FOL.m s function [ y ] = TF_FOL(s,para) y = 1./(ones(size(s))+para(1).*s); s G( s) = 1 1+ τs Matlab

25 FOL_boad_ma in.m function main tau = 1; % j=sqrt(-1);% para = [tau];% omega = 10.^(-2:0.1:2); % G = TF_FOL(j.*omega,para);% gain = 20.*log10(abs(G)); % kaku = angle(g); % subplot(2,1,1); %2 semilogx(omega,gain); %x xlabel( omega ); ylabel( db );% subplot(2,1,2); %2 semilogx(omega,kaku); %x xlabel( omega ); ylabel( rad ); % end

26 MATLAB Runge-Kutta

27 d 2 x( t) dt 2 dx( t) + 2ζω 2 n dt + ω 2 x( t) = ω 2 u( t) n n xx (0)=0 L{ u( t) } = U( s) L{ x( t) } = X( s) ( ) = X ( s ) U( s) = ω 2 n G s s 2 + 2ζω n s + ω n 2

28 U( s) = 1 s X( s) = ω n 2 s 2 + 2ζω n s + ω n 2 1 s 3 ζ <1 ζ =1 ζ >1

29 ζ <1 X s ( ) = 1 s ω 2 n 2 1 β β + jα ( ) α = ω n ζ 1 s + α jβ ω 2 n 2 1 β β jα ( ) β = ω n 1 ζ 2 1 s + α + jβ x( t) =1 1 1 ζ exp ( ω nζt)sin 1 ζ 2 ω 2 n t + φ ( ) φ = tan 1 1 ζ 2 ζ ζ =1 X( s) = 1 s 1 s + ω n ω n ( ) 2 s + ω n x( t) =1 ( 1+ ω t)exp( ω t) n n

30 ζ >1 X( s) = 1 s + s s 1 s 2 s s 1 s 1 1 s 2 s 1 s s 2 s 1 = ζω n + ω n ζ 2 1 s 2 = ζω n ω n ζ 2 1 x( t) =1+ s 2 exp( s 1 t) + s 1 s 2 s 1 exp( s 2 t) s 2 s 1

32 MATLAB Runge-Kutta

33 2 d 2 x t ( ) dt 2 ( ) 2 dx t + 2ζω n dt + ω n 2 x t ( ) = ω n 2 u t ( ) dx t ( ) dt = y t ( ) d 2 x t ( ) dt 2 = dy ( t ) dt = 2ζω n 2 y t ( ) ω n 2 x t ( ) + ω n 2 u t ( ) d dt x = y ω n 2ζω n x y + 0 ω n 2 u t ( ) x = [ 1 0] x y

34 Runge-Kutta RK2.m SOE.m SOE_main.m

35 (SOE.m) d dt x = y y ω 2 n x 2ζω 2 n y + ω 2 n u t ( )

36 SOE_main.m function main xi = 0.2; % wn = 2; % In = 1;% dt = 0.1;% T = 15;% xinit =[0;0];% para = [xi,wn,in];% %@SOE [t1,xy] = RK2(@SOE,[0,T],xinit,dt,para); %2 Runge-Kutta x = xy(1,:); xmin = min(x)-0.1; xmax = max(x)+0.1; y = xy(2,:); ymin = min(y)-0.1; ymax = max(y)+0.1; subplot(2,1,1); %2 plot(t1,x);% axis([0 T xmin xmax]); xlabel('t'); ylabel('x'); subplot(2,1,2); %2 plot(x,y);% axis([xmin xmax ymin ymax]); xlabel('x'); ylabel('dx/dt'); end

37 MATLAB Runge-Kutta

38 Runge-Kutta RK2.m SOE_FR.m SOE_FR_main.m

39 (SOE.m) function xy = SOE(t,y,para) xy = [y(2); % -para(2).^2.*y(1)-2.*para(1).*para(2).*y(2)+para(2).^2.*para(3).*sin(para(4).*t)]; end d dt x y = y ω 2 n x 2ζω 2 n y + ω 2 n Asin ωt ( )

40 SOE_FR_mai n.m xi = 0.5; % wn = 1; % In = 1; % omega = 1.5; % dt = 0.01; % T = 60; % xinit =[0;0]; % para = [xi,wn,in,omega]; % %@SOE_FR [t1,xy] = RK2(@SOE_FR,[0,T],xinit,dt,para); %2 Runge- Kutta x = xy(1,:); xmin = min(x)-0.1; xmax = max(x)+0.1; y = xy(2,:); ymin = min(y)-0.1; ymax = max(y)+0.1; subplot(2,1,1); %2 plot(t1,x); % axis([0 T xmin xmax]); xlabel('t'); ylabel('x'); subplot(2,1,2); %2 plot(x,y); % axis([xmin xmax ymin ymax]); xlabel('x'); ylabel('dx/dt');

41 d 2 x t ( ) dt 2 ( ) 2 dx t + 2ζω n dt + ω n 2 x t ( ) = ω n 2 sin ωt ( ) ( ) x( t) = G( jω) sin ωt + tan 1 G( jω) G( jω) = = ω n 2 ω jζω ω + ω 2 n n 1 ω 2 ω jζ ω ω +1 n n

42 ω ω n

43 TF_SOE.m FOL_boad_main.m

44 TF_SOE.m s function [ y ] = TF_SOE(s,para) y = para(2).^2./(s.^2 + 2.*para(1).*para(2).*s + para(2).^2.*ones(size(s))); end s G( s) = ω n 2 s 2 + 2ζω n s + ω n 2 Matlab

45 SOE_boad_ma in.m function main xi = 0.05;% wn = 1;% j=sqrt(-1);% para = [xi,wn]; % omega = 10.^(-1:0.05:1); % G = TF_SOE(j.*omega,para); % gain = 20.*log10(abs(G)); % kaku = angle(g); % subplot(2,1,1); %2 semilogx(omega,gain); %x xlabel('omega'); ylabel('db'); % subplot(2,1,2); %2 semilogx(omega,kaku); %x xlabel('omega'); ylabel('rad'); % end

46 RK2.m 4 Runge-Kutta 4 Runge-Kutta 2

47 4 Runge-Kutta x n +1 = x n + Δt t = nδt [ ] 6 F 1 + 2F 2 + 2F 3 + F 4 F = f ( x,t) 1 n F = f x + Δt 2 n 2 F Δt,t F = f x + Δt 3 n 2 F Δt,t F = f ( x + ΔtF,t + Δt) 4 n 3

x () g(x) = f(t) dt f(x), F (x) 3x () g(x) g (x) f(x), F (x) (3) h(x) = x 3x tf(t) dt.9 = {(x, y) ; x, y, x + y } f(x, y) = xy( x y). h (x) f(x), F (x

x () g(x) = f(t) dt f(x), F (x) 3x () g(x) g (x) f(x), F (x) (3) h(x) = x 3x tf(t) dt.9 = {(x, y) ; x, y, x + y } f(x, y) = xy( x y). h (x) f(x), F (x [ ] IC. f(x) = e x () f(x) f (x) () lim f(x) lim f(x) x + x (3) lim f(x) lim f(x) x + x (4) y = f(x) ( ) ( s46). < a < () a () lim a log xdx a log xdx ( ) n (3) lim log k log n n n k=.3 z = log(x + y ),