% This file allows to plot some nonlinear systems: % Lorenz, Duffing, Chua circuit, Colpitts and Rayleigh % oscillators, Van der Pol oscillator (for µ=0.2, µ=1, µ=10), % and Henon map. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % [1] Plot Lorenz dynamic by using ode45 command: % % Lorenz's equations: % % dx/dt = sigma*(y-x) % dy/dt = x*(rho-z)-y % dz/dt = x*y-beta*z % % with sigma=10, rho=28., beta=8./3, % Initial conditions x0=[10 10 10]. disp('[1] Lorenz ') disp(' ') disp(' ... Please Wait ...') disp(' ') disp('Then Press any Key to Continue ...') disp(' ') disp(' ') global sigma r b sigma=10.; r=28.; b=8./3.; x0=[10 10 10]; t0=0; tf=100; [t, x] = ode45(@lorenzeqs, [t0, tf], x0); figure; plot3(x(:,1),x(:,2),x(:,3)); title('Lorenz'); xlabel('x'); ylabel('y'); zlabel('z'); axis([-10 10 -30 30 0 70]); box; pause clear all %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % [2] Plot Duffing dynamic in the Phase space by using ode45 command % % Duffing's equations (general form): % % du/dt = v % dv/dt = -(omega_0)^2*u-beta*u^3-delta*v+gamma*cos(omega*t+phi) % % where u is the displacement, v is the velocity, and the % parameters are constants. disp('[2] Duffing') disp(' ') disp(' ... Please Wait ...') disp(' ') disp('(To Stop the Simulation, Press "Ctrl+c" in the Command Window.)') disp(' ') disp('Press any Key to Continue, ....') disp(' ') [t,z]=ode45(@duffing,[0,300],[3.01,4.01]); plot(z(:,1),z(:,2)); title('Duffing, phase space'); pause clear all %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % [3] Plot Chua circuit dynamic by using ode45 command % % Chua's autonomous circuit % % dx/dt = alpha*(-x+y-H(x)) % dy/dt = x-y+z % dz/dt = -beta*y % % where H(x) = bx+0.5*(a-b)*(|x+1|-|x-1|), % with alpha arbitrary, alpha is selected here equal to 10, % beta = 15, a = -1.3, b = -0.7 disp('[3] Chua circuit') disp(' ') disp(' ... Please Wait ...') disp(' ') disp('(To Stop the Simulation, Press "Ctrl+c" in the Command Window.)') disp(' ') disp('Press any Key to Continue, ....') disp(' ') y0=[1 0.2 0]; [t,z]=ode45(@Chua,[0 300],y0); plot(z(:,1),z(:,2)); title('Chua Circuit'); pause clear all %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % [4] Plot Colpitts Oscillator dynamic by using ode45 command % % Colpitts oscillator system: % % dx/dt = z-beta*f(y) % dy/dt = alpha*(-r*(y+gamma)-z-f(y)) % dz/dt = delta*(-x+y+epsilon)-rho*z % % where f(y)=0 if y<=1, and y-1 otherwise. % with alpha=1, beta = 200, gamma = -20/3, delta = 5.5, % with r arbitrary, rho = 2,epsilon = 20/3. disp('[4] Colpitts') disp(' ') disp(' ... Please Wait ...') disp(' ') disp('Then Press any Key to Continue, ....') disp(' ') y0=[0.1 0 0]; s=30;Tm=[0:1/s:200]; [t,y]=ode45(@Colpitts,Tm,y0); plot3(y(:,1),y(:,2),y(:,3)); xlabel('Y_1'); ylabel('Y_2'); zlabel('Y_3'); title('Colpitts Oscillator'); axis([0 8 -2.5 2 -2 8]) view(70,25); box; pause clear all %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % [5] Rayleigh oscillator with physical parameters % % c*dv/dt = i+id(-v) % b*di/dt = -v-u1-u2*sin(omega*t) % % where u1=1, u2=878, omega=60.3e3, a1=1.542, a2=0.5, % b=10, c=2.377 and id(v)=-a1*v+a2*(-v)^3. disp('[5] Rayleigh oscillator') disp(' ') disp(' ... Please Wait ...') disp(' ') disp('(To Stop the Simulation, Press "Ctrl+c" in the Command Window.)') disp(' ') disp('Press any Key to Continue, ....') disp(' ') y0=[0.1 0]; q=10e5; vtm=[0:1/q:20e-3]; [t,y]=ode45(@Rayleigh,vtm,y0); plot(y(:,1),y(:,2)); xlabel('y_1'); ylabel('y_2'); title('Rayleigh Oscillator'); pause clear all %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % [6] The van der Pol equation (y_1)"-µ*(1-(x_1)^2)*y'+(y_1)=0 is equivalent % to a system of coupled first-order differential equations % (y_1)'=(y_2) % (y_2)'=µ*[1-(y_1)^2]*(y_2)-(y_1). % Evaluate the van der Pol Ode for mu = 0.2. % The ode command is as follows: disp('[6] van der Pol oscillator for µ=0.2') disp(' ') disp(' ... Please Wait ...') disp(' ') disp('(To Stop the Simulation, Press "Ctrl+c" in the Command Window.)') disp(' ') disp('Press any Key to Continue ...') disp(' ') x0=[0.1 0.1]; [t,y] = ode45(@vdpeq1,[0 100],x0); plot(t,y(:,1),'-',t,y(:,2),'-.'); xlabel('time t'); ylabel('solution y_1 and y_2'); title('van der Pol with \mu=0.2'); pause disp(' ') disp('[ . 6 bis] ') disp(' ') disp('Press any Key to Continue ...') disp(' ') plot(y(:,1),y(:,2)); xlabel('y_1'); ylabel('y_2'); title('van der Pol with \mu=0.2'); pause clear all %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % [7] The van der Pol equation (y_1)"-µ*(1-(x_1)^2)*y'+(y_1)=0 is equivalent % to a system of coupled first-order differential equations % (y_1)'=(y_2) % (y_2)'=µ*[1-(y_1)^2]*(y_2)-(y_1). % Evaluate the van der Pol Ode for mu = 1. % The ode command is as follows: disp('[7] van der Pol oscillator for µ=1') disp(' ') disp('Press any Key to Continue ....') disp(' ') x0=[2 0]; [t,y] = ode45(@vdpeq2,[0 20],x0); plot(t,y(:,1),'-',t,y(:,2),'-.'); xlabel('time t'); ylabel('solution y_1 and y_2'); title('van der Pol with \mu=1'); pause disp(' ') disp('[ . 7 bis] ') disp(' ') disp('Press any Key to Continue ....') disp(' ') plot(y(:,1),y(:,2)); xlabel('y_1'); ylabel('y_2'); title('van der Pol with \mu=1') pause clear all %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % [8] The van der Pol equation (y_1)"-µ*(1-(x_1)^2)*y'+(y_1)=0 is equivalent % to a system of coupled first-order differential equations % (y_1)'=(y_2) % (y_2)'=µ*[1-(y_1)^2]*(y_2)-(y_1). % Evaluate the van der Pol Ode for mu = 10. % The ode command is as follows: disp(' ') disp('[8] van der Pol oscillator for µ=10') disp(' ') disp('Press any Key to Continue, ....') disp(' ') x0=[2 0]; [t,y] = ode15s(@vdpeq3,[0 100],x0); plot(t,y(:,1),'-'); xlabel('time t'); ylabel('solution y_1'); title('van der Pol with \mu=10') pause disp(' ') disp('[ . 8 bis] ') disp(' ') disp('Press Any Key to Continue ...') disp(' ') plot(y(:,1),y(:,2)); xlabel('y_1'); ylabel('y_2'); title('van der Pol with \mu=10') pause clear all %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % [9] Plot the Henon map by using the function henonmap(iterations); % (The number of iterations has to be selected between 100 and 2000). disp('[9] Henon map') disp(' ') disp(' ... Please Wait Few Seconds ...') disp(' ') disp('(To Stop the Simulation, Press "Ctrl+c" in the Command Window.)') disp(' ') henonmap(1500); % "Complex and Chaotic Nonlinear Dynamics. % Advances in Economics and Finance, % Mathematics and Statistics" % T.Vialar, Springer 2009 % Copyright(c).