sol_ex1_ss18

ex1_sol.pdf PDF document, 423 KB

code_ex11

Exercise1_1.m Objective-C source code, 2 KB

File contents

function Exercise1_1

 % Solution to Exercise 1.1
 
 % Initialization
 starttime = 0;                          % Starttime
 endtime   = 5;                          % Endtime
 deltaT    = 0.1;                        % Timesteps
 tsim      = (starttime:deltaT:endtime); % Timevector
 
 x0      = 5;                            % Inital state value
 lambda = -1;                            % Additional Paramter
 
 % Inital States for Euler, Improved Euler, Runge-Kutta, analytical method
 x_Euler1    = x0;
 x_impEuler1 = x0;
 x_4RK1      = x0;
 x_an1       = x0;
 
 % Inital States for plot
 xs_1(1) = x0; % Euler
 xs_2(1) = x0; % Improved Euler
 xs_3(1) = x0; % Runge-Kutta
 xs_4(1) = x0; % Analytical
 
 % Numerical Simulation
 for i = 1:1:length(tsim)-1;    % Loop timevector length-1 because t0
                                % is allready initialized
   
     % Analytical solution
     x_an2 = x_an1*exp(lambda*i*deltaT);
    
     % Euler's method
     x_Euler2 = x_Euler1 + lambda*x_Euler1*deltaT;     
     
     % Improved Euler method
     x_trial     = x_impEuler1 + lambda*x_impEuler1*deltaT;
     x_impEuler2 = x_impEuler1 + 1/2*(lambda*x_impEuler1 + lambda*x_trial)*deltaT;
     
     % Fourth-order Runge-Kutta method.
     k1 = lambda*x_4RK1*deltaT;
     k2 = lambda*(x_4RK1+0.5*k1)*deltaT;
     k3 = lambda*(x_4RK1+0.5*k2)*deltaT;
     k4 = lambda*(x_4RK1+0.5*k3)*deltaT;
     
     x_4RK2 = x_4RK1 + (1/6 * (k1 + 2*k2 + 2*k3 + k4));
     
     % Prepare next iteration
     x_Euler1    = x_Euler2;
     x_impEuler1 = x_impEuler2;
     x_4RK1      = x_4RK2;
     
     % Store for plot
     xs_1(i+1) = x_Euler1;
     xs_2(i+1) = x_impEuler1;
     xs_3(i+1) = x_4RK1;
     xs_4(i+1) = x_an2;
     
 end  % End simulation loop

 % Alternative using predefined ode45-solver
 % Call ode45-solver (see "doc ode45" for more information)
 [t1,y1] = ode45(@(t,x)exfun(t,x,lambda),tsim,x0);
 
 % System defined by ODEs  [outputs] = functionname(inputs)
 function dx = exfun(t,x,lambda)
 
     dx = lambda*x;
     
 end

 % Plot results of Euler, improved Euler, Runge Kutta, analytical, ode45
 figure(1)
 plot(tsim,xs_1,tsim,xs_2,tsim,xs_3,tsim,xs_4,t1,y1,'LineWidth',2)
 grid on
 legend('Euler','improved Euler','4th-Runge-Kutta','Analytical','ode45')

end

code_ex14

Exercise1_4.m Objective-C source code, 2 KB

File contents

function Exercise1_4

 % Solution to Exercise 1.4
 
 % Initialization
 starttime = 0;                          % Starttime
 endtime   = 10;                         % Endtime
 deltaT    = 0.1;                        % Timesteps
 tsim      = (starttime:deltaT:endtime); % Timevector

 % Inital state value
 x0 = [1, 3];                 
 
 % System matrix a)
 Aa = [-4, 2; -3,  1];                   % Additional Paramter
 % b)
 Ab = [-1, 1; -1, -1];
 % b2)
 Ab2 = Ab;
 % c)
 Ac = [-2, 1; -1, -4];
 % d)
 Ad = [ 2, 0;  5, -1];
 
 % Store initial values for plot
 xs_1a(1) = x0(1);               % x1
 xs_2a(1) = x0(2);               % x2
 xs_1b(1) = x0(1);
 xs_2b(1) = x0(2);
 xs_1c(1) = x0(1);
 xs_2c(1) = x0(2);
 xs_1d(1) = x0(1);
 xs_2d(1) = x0(2);
 
 % Simulation loop
 for i = 1:1:length(tsim)-1;    % Loop timevector length-1 because t0
                                % is allready initialized
   
     % Analytical solution
     % System a)
     x_an12a = (3*exp(-2*deltaT*i)-2*exp(-deltaT*i))*x0(1)...
               +(-2*exp(-2*deltaT*i)+2*exp(-deltaT*i))*x0(2);
           
     x_an22a = (3*exp(-2*deltaT*i)-3*exp(-deltaT*i))*x0(1)...
               +(-2*exp(-2*deltaT*i)+3*exp(-deltaT*i))*x0(2);
     
     % System b)
     x_an12b = exp(-deltaT*i)*(cos(deltaT*i)*x0(1)+sin(deltaT*i)*x0(2));
     x_an22b = exp(-deltaT*i)*(-sin(deltaT*i)*x0(1)+cos(deltaT*i)*x0(2));
     
     % System c)
     x_an12c = (1+deltaT*i)*exp(-3*deltaT*i)*x0(1)+(deltaT*i)*exp(-3*deltaT*i)*x0(2);
     x_an22c = (-deltaT*i)*exp(-3*deltaT*i)*x0(1)+(1-deltaT*i)*exp(-3*deltaT*i)*x0(2);
     
     % System d)
     x_an12d = exp(2*deltaT*i)*x0(1);
     x_an22d = (5/3)*(exp(2*deltaT*i)-exp(-deltaT*i))*x0(1) + exp(-deltaT*i)*x0(2);
           
     % Store for plot a)b)c)d)
     xs_1a(i+1) = x_an12a;
     xs_2a(i+1) = x_an22a;
     xs_1b(i+1) = x_an12b;
     xs_2b(i+1) = x_an22b;
     xs_1c(i+1) = x_an12c;
     xs_2c(i+1) = x_an22c;
     xs_1d(i+1) = x_an12d;
     xs_2d(i+1) = x_an22d;
      
 end
 
 % Call ode45-solver (see "doc ode45" for more information)
 % System Aa
 [ta,ya] = ode45(@(t,x)exfun(t,x,Aa),tsim,x0);
 % System Ab
 [tb,yb] = ode45(@(t,x)exfun(t,x,Ab),tsim,x0);
 % System Ac
 [tc,yc] = ode45(@(t,x)exfun(t,x,Ac),tsim,x0);
 % System Ad
 [td,yd] = ode45(@(t,x)exfun(t,x,Ad),tsim,x0);

 
 % System defined by ODEs  
 function dx = exfun(t,x,A)
 
     dx = A*x;
     
 end

% Plot results
figure(1)
subplot 221
plot(ta,ya,tsim,xs_1a,tsim,xs_2a,'LineWidth',2)
title('a)')
grid on
legend('x_1_o_d_e_4_5','x_2_o_d_e_4_5','x_1_a_n_a_l_y_t','x_2_a_n_a_l_y_t')
subplot 222
plot(tb,yb,tsim,xs_1b,tsim,xs_2b,'LineWidth',2)
title('b)')
grid on
legend('x_1_o_d_e_4_5','x_2_o_d_e_4_5','x_1_a_n_a_l_y_t','x_2_a_n_a_l_y_t')
subplot 223
plot(tc,yc,tsim,xs_1c,tsim,xs_2c,'LineWidth',2)
title('c)')
legend('x_1_o_d_e_4_5','x_2_o_d_e_4_5','x_1_a_n_a_l_y_t','x_2_a_n_a_l_y_t')
grid on
subplot 224
plot(td,yd,tsim,xs_1d,tsim,xs_2d,'LineWidth',2)
title('d)')
legend('x_1_o_d_e_4_5','x_2_o_d_e_4_5','x_1_a_n_a_l_y_t','x_2_a_n_a_l_y_t')
grid on

end

solEx32

Ex32_ss2018.m Objective-C source code, 2 KB

File contents

function Ex22(x0)
	% octave code for solving simultaneously the time-varying ode
	% and the matrix differential equation for the fundamental matrix
	% as well as the monodromy matrix, Floquet multipliers and Lyapunov
	% exponents.
	% Call with Ex22(x0)
	% where x0 is the initial condition for the simulation.
	% A. Schaum, 2018


	% Simulation time vector
	T = 2*pi % 1 period of the matrix A
	t = linspace(0,T,200);

	% Initial condition for the fundamental matrix in vector form,
	% chosen as identity for the purpose of calculating the monodromy
	% matrix and the Floquet exponents
	p0 = [1;0;0;1];
    
	% Solving the differential equations
	[ts, sol] = ode45(@dyn,t,[x0;p0]);

	% Determination of the monodromy matrix if the simulation time T
	% has been chosen as one period!
	M = [sol(end,3), sol(end,4);
	     sol(end,5), sol(end,6)];
	[v,lamb] = eig(M);
	mu=[lamb(1,1);lamb(2,2)];
	
	% Determining the logarithm of M (Notice, that M is diagonalizable
	% with diagonal form M1..)

	M1 = inv(v)*M*v;            % Diagonalization of M
	B=zeros(2,2);               % Initialization of B = 1/T*log(M)
	B(1,1) = 1/T*log(M1(1,1));  
	B(2,2) = 1/T*log(M1(2,2));
	sigma = real(eig(B));       % Lyapunov exponents

	fprintf('Monodromy matrix:\n');
	M
	fprintf('Floquet multipliers:\n');
	mu
	fprintf('log(M):\n');
	B
	fprintf('Lyapunov exponents:\n');
	sigma

	% Graphical output
	% Time responses of the state variables
	figure(1)
	subplot(2,1,1)
	plot(ts,sol(:,1))
	ylabel('x_1')
	grid on
	subplot(2,1,2)
	plot(ts,sol(:,2))
	ylabel('x_2')
	xlabel('Time t')
	grid on
    	% Phase plane plot
	figure(2)
	plot(sol(:,1),sol(:,2))
	grid on
	hold on

	

end

function dxdt = dyn(t,x)
	u=0;% place holder for control problem..

	% Changing back to the matrix notation for P
	P=[x(3),x(4); 
	   x(5),x(6)];
	% the system dynamics matrix
	A = [-1+cos(t), 0; 1, -.5+sin(t)];
	% Determining the right-hand side of the ode 
	dxdt(1:2,1) = A*x(1:2)+[1;0]*u;	

	% and the matrix differential equation:
	dPdt = A*P;

	% completing the right-hand-side in vector notation
	dxdt(3,1) = dPdt(1,1);
	dxdt(4,1) = dPdt(1,2);
	dxdt(5,1) = dPdt(2,1);
	dxdt(6,1) = dPdt(2,2);

end

E10_1.m

Ex10_1.m Objective-C source code, 1 KB

File contents

function Ex10_1(N,x0)

xv(1,1) = x0;
xv(2,1) = 0;
x=x0;
for i = 2 : N
  xv(1,2*(i-1))= x;
  xv(2,2*(i-1))= f(x);
  xv(1,2*(i-1)+1)= f(x);
  xv(2,2*(i-1)+1)= f(x);
  x = f(x);
end
xn=linspace(-1.2,1.2,100);
figure;
plot(xv(1,:), xv(2,:),'r'); hold on;
plot(xn,f(xn));hold on
plot(xn,xn,'k');
%axis([0,1.2,0,1.2])
grid on

figure
grid on; hold on;
plot(xv(1,1:2:end),'*',xv(1,1:2:end),'m')

end

function r=f(x)
  r= -1/2*x+x.^3;
end

Ex10_4.m

Ex10_4.m Objective-C source code, 1 KB

File contents

function X=Ex10_4(N,xc0)
xv(1,1) = xc0(1);
xv(2,1) = 0;
x=xc0(1);
Xa = xc0;
X(:,1) = Xa;
for i = 2 : N
  %Cobweb plot solution
  xv(1,2*(i-1))= x;
  xv(2,2*(i-1))= f_c(x);
  xv(1,2*(i-1)+1)= f_c(x);
  xv(2,2*(i-1)+1)= f_c(x);
  x = f_c(x);
  %State-space solution
  X(:,i) = f(Xa);
  Xa=X(:,i);
end

xca=linspace(-0.2,0.2,100);
% Cobweb on the center manifold
figure(1);
plot(xv(1,:),xv(2,:),'r');hold on;
plot(xca,f_c(xca));%hold on
plot(xca,xca,'k');
grid on
% Time response on the center manifold
figure(2)
plot(xv(1,1:2:end),'*');hold on;
plot(xv(1,1:2:end),'m')
grid on; hold on;
%State space plot
figure(3)
x2n = linspace(-0.5,0.5,100);
x1n = 2*x2n.^2; %approx. center manifold
plot(x1n,x2n,X(1,1:10:end),X(2,1:10:end),'*')
hold on
end

function out = f(Xa)
  x1 = Xa(1);
  x2 = Xa(2);
  
  out = [(1/2-x2)*x1+(x1+x2)^2;
              -x2+x1^3];
end

function out = f_c(x2)
  x1 = 2*x2;
  out = -1*x2+x1.^3;
end
Research