preface

The following are some boutique columns I have prepared for you. My favorite partners can subscribe by themselves. Your support is the driving force for me to update!

MATLAB-30 days take you from introduction to mastery

Advanced course for in-depth understanding of MATLAB (with source code)

tableau visual data analysis advanced tutorial

python fast learning practical application series courses

Nonlinear phase diagram

We use several examples to introduce the drawing method of phase diagram of nonlinear system.

The examples given here are all examples of autonomous system equations, that is, the result of the equation is independent of the initial value of t0 (time invariant system), and does not contain quantities related to t such as external periodic driving force. Because to describe an autonomous system, we only need to know the derivatives of various variables in the space of the system, and then form the phase space. However, each state of the time-varying system will change with time, so the static phase diagram can not be used for qualitative analysis.

Therefore, for autonomous second-order systems, the second-order phase plane can completely describe the motion state of the system, whether linear or nonlinear. The usual steps can be divided into two steps: (1) calculate the dy and ddy derivatives of each point, (2) draw the streamline diagram corresponding to the vector field according to the vector obtained at each point.

Take an example given in the book nonlinear systems as an example. The formula of the second-order nonlinear equation is as follows:

The phase diagram of the nonlinear system can be obtained by drawing the system derivative of each point in space and drawing the streamline.

It can be seen that the phase plane of a nonlinear system may have more than one equilibrium point. However, the neighborhood near these equilibrium points can also be linearized, and its behavior can be explained by the phase diagram of the linear system in the previous section. For example, the figure above has two stable nodes and a saddle point.

Next, we introduce a classical phase plane pattern: limit cycle, which only appears under nonlinear conditions.

Taking the classical Van der Pol equation as an example, the form of this equation is as follows:

hinder ε Is a constant, ε The larger the, the greater the nonlinearity of the equation. take ε= 0.3, draw the phase space diagram:

There is only one equilibrium point in the whole space, that is, the spiral point divergent in the limit cycle. In order to intuitively understand the characteristics of limit cycle, different initial conditions are selected to solve the equation in time domain, and the following results are obtained:

At the point outside the limit cycle, it initially presents a convergent waveform, but then the amplitude remains constant. At the point in the limit cycle, it initially presents a divergent waveform, but then the amplitude remains stable. Both of them are attracted by the limit cycle in space and carry out circle drawing movement again and again.

If some coefficients in the equation change (such as those in the above equation) ε Changing from negative to positive) will change the whole phase space. This change will change the properties of the equilibrium point and change the properties of the whole system. We call the system bifurcation at this position. For example, the Van der Pol equation ε From negative to positive, the equilibrium point changes from stable to divergent, resulting in the stable position in space from a point to a limit cycle of motion.

If the vibration signal can be observed in the experiment, the phase plane can also be drawn to observe the characteristics of the signal. The phase plan of common periodic signals in the measurement process is given in the following three diagrams:

The first figure is a typical linear vibration diagram, the waveform is a sinusoidal curve, and the phase plane is a circle.

The second picture shows a typical limit cycle vibration. One characteristic of this kind of vibration is that even if disturbed, the amplitude can always remain at a stable value. If this phenomenon is found in practical experiments and there is no external periodic driving force (the equation is an autonomous system), it can be basically determined as limit cycle motion.

The third picture shows a typical high-dimensional nonlinearity. Because the flow lines in the phase plane do not cross. This cross curve is the projection of high-dimensional space on a two-dimensional plane. The figure shows the simulation of period doubling phenomenon in high-dimensional nonlinearity. This will be introduced in later articles.

Of course, due to the noise in the actual test, the measurement results will not be so good-looking. For example, if the final result is chaotic and irregular, chaos must appear in theory, but in fact, it is likely that the noise is too large.

The matlab code used in the drawing is attached below:

%1 Two dimensional phase space%Nonlinear clearclcclose all
%1 Nonlinear systems with multiple equilibria%reference resources Nonlinear system (third edition of Chinese translation) Khalil P32[y,dy]=meshgrid(-0.5:0.02:1.5,-0.5:0.02:1.5);%Initialize grid u=zeros(size(y));v=u;for k=1:numel(y)    %Calculate the direction of each point on the grid    F=Fdydx(0,[y(k);dy(k)],1);    u(k)=F(1);    v(k)=F(2);endfigure()streamslice(y,dy,u,v,4)xlabel('y')ylabel('dy')box on
%2 Limit cycle[y,dy]=meshgrid(-4:0.1:4,-4:0.1:4);%Initialize grid u=zeros(size(y));v=u;for k=1:numel(y)    %Calculate the direction of each point on the grid    F=Fdydx(0,[y(k);dy(k)],2);    u(k)=F(1);    v(k)=F(2);endfigure()streamslice(y,dy,u,v,3)xlabel('y')ylabel('dy')box onxlim([-4,4]);ylim([-4,4])
figure()streamslice(y,dy,u,v,3)xlabel('y')ylabel('dy')box onhold ont2=0:0.1:50;[y2,~]=ODE_RK4_hyh(t2,t2(2)-t2(1),[-4;4],2);plot(y2(1,:),y2(2,:),'r','linewidth',2)[y3,~]=ODE_RK4_hyh(t2,t2(2)-t2(1),[-0.1;0.1],2);plot(y3(1,:),y3(2,:),'g','linewidth',2)hold offxlim([-4,4]);ylim([-4,4])figure()plot(t2,y2(1,:),'r','linewidth',2);ylim([-4,4])figure()plot(t2,y3(1,:),'g','linewidth',2);ylim([-4,4])
%3 How different time domain signals look on the phase diagram figure()t1=0:0.01:1;y1=sin(10*pi*t1);dy1=diffhyh(y1,2);subplot(2,3,1)plot(t1,y1);ylim([-1.5,1.5])xlabel('t');ylabel('y')subplot(2,3,4)plot(y1,dy1);xlim([-1.5,1.5])xlabel('y');ylabel('dy')%Nonlinear limit cycle
t2=0:0.01:18;[y,Output]=ODE_RK4_hyh(t2,t2(2)-t2(1),[-1;1.161],2);y2=y(1,:);dy2=y(2,:);subplot(2,3,2)plot(t2,y2);ylim([-2.5,2.5]);xlim([0,18])xlabel('t');ylabel('y')subplot(2,3,5)plot(y2,dy2);xlim([-2.5,2.5])xlabel('y');ylabel('dy')
%Chaos (simulated signal, no actual system)%In the two-dimensional system, the phase plane can completely describe the system, and there can be no intersection.%If there is trajectory crossing, it means that the dimension of the system is greater than 2, and the graph displayed by the phase plane is only the projection of the high-dimensional system on the 2 dimension. t3=0:0.01:2;Gas=@(t,a,b,c) c*exp(-b*(t-a).^2);y3=Gas(t3,0,120,1)+Gas(t3,1,120,1)+Gas(t3,2,120,1);y3=y3+Gas(t3,0.3,90,0.5)+Gas(t3,1.3,90,0.5);y3=y3-Gas(t3,0.6,300,0.3)-Gas(t3,1.6,300,0.3);dy3=diffhyh(y3,2);subplot(2,3,3)plot(t3,y3);ylim([-0.5,1.5])xlabel('t');ylabel('y')subplot(2,3,6)plot(y3,dy3);xlim([-0.5,1.5]);ylim([-0.15,0.15])xlabel('y');ylabel('dy')
function [F,Output]=Fdydx(x,y,Input)%Form Y'=F(x,Y)For the equation, see the relevant knowledge of numerical analysis and solving constant coefficient differential equations%The higher order is represented by a column vector, F=[dy(1);dy(2)];y(1)Is a function, y(2)Is the derivative of the function switch Input    %Example 1    case 1        p=[83.72,-226.31,229.62,-103.79,17.76,0];        dy(1)=0.5*(-polyval(p,y(1))+y(2));        dy(2)=0.2*(-y(1)-1.5*y(2)+1.2);        F=[dy(1);dy(2)];    %Example 2    case 2        dy(1)=y(2);        dy(2)=-y(1)+0.3*(1-y(1)^2)*y(2);%0.3        F=[dy(1);dy(2)];endOutput=[];end
function [y,Output]=ODE_RK4_hyh(x,h,y0,Input)%4 rank RK method%h Algorithm with constant interval y=zeros(size(y0,1),size(x,2));y(:,1)=y0;for ii=1:length(x)-1    yn=y(:,ii);    xn=x(ii);    [K1,~]=Fdydx(xn    ,yn       ,Input);    [K2,~]=Fdydx(xn+h/2,yn+h/2*K1,Input);    [K3,~]=Fdydx(xn+h/2,yn+h/2*K2,Input);    [K4,~]=Fdydx(xn+h  ,yn+h*K3  ,Input);    y(:,ii+1)=yn+h/6*(K1+2*K2+2*K3+K4);endOutput=[];end
function Fdiff=diffhyh(F,dim)%The second-order central difference is adopted, and the first-order forward or backward difference is adopted for the edge (there is no upwind difference at the edge, and the accuracy requirement is not high)%dim，The dimensional direction of the difference, dim=1 Corresponding to the difference in the downward direction of the matrix, dim=2 Corresponding right diffF1=diff(F,1,dim);if dim==1%down    Fdiff=([zeros(1,size(F,2));diffF1]+[diffF1;zeros(1,size(F,2))])/2;    Fdiff(1,:)=diffF1(1,:);    Fdiff(end,:)=diffF1(end,:);elseif dim==2%towards the right    Fdiff=([zeros(size(F,1),1),diffF1]+[diffF1,zeros(size(F,1),1)])/2;    Fdiff(:,1)=diffF1(:,1);    Fdiff(:,end)=diffF1(:,end);endend

Chaotic system

If there are intersecting trajectories in the two-dimensional phase plane, it indicates that the dimension of the system is likely to be greater than two-dimensional.

The following is a demonstration of several classic systems. This chapter does not involve too many knowledge points, mainly display. Three classical nonlinear chaotic systems are introduced.

• one   Lorenz system

Lorenz system is a nonlinear system discovered and proposed by Lorenz, a meteorologist. It is also the beginning of chaos. When simulating atmospheric flow, Lorentz found that a small initial error will lead to great changes in the system in the future. This thought dealt a heavy blow to the determinists in the field of physics in the 1960s. Lorenz also summarized this uncertainty as "Butterfly Effect".

This system can be written as:

Generally, the system a=10, b=8/3. Change the r value to observe the different appearance of the system. The following figure shows the two-dimensional phase trajectory diagram of xy plane corresponding to different r values:

Where r=20, the corresponding system converges to a fixed point. r=28 corresponds to chaos. r=99.36 corresponds to period doubling.

Because the two-dimensional system will not cross on the phase plane, chaos, period doubling and other phenomena only appear in three-dimensional or higher dimensions. As mentioned above, for chaos, three is a magical number.

Their trajectories in three-dimensional space are:

The graph in the middle is the classical Lorentz attractor graph.

• 2 Rossler system

Rossler system is a nonlinear system proposed by R ö ssler himself in the 1970s. Compared with the previous Lorenz system, Rossler system is simpler, but it still has complex nonlinear behavior.

It can write:

In the following figure, a=0.1 and b=0.1 are plotted, and different c plots are changed.

Cycle 2 refers to a cycle for every two waveforms, and the system turns 2 times to return to the same point. With the increase of c, the system changes from period 1 to period 2, then suddenly increases large period 4, and then increases to period 8 or even higher at a faster and faster speed. Finally, the dense period seems to never cycle and become chaos.

This phenomenon that the period becomes chaotic gradually is called period doubling bifurcation. It is a common way for a system to change from order to disorder and chaos.

• 3 duffing equation

duffing equation is also based on   A nonlinear equation named by Georg Duffing. It is based on the equation of forced vibration simple pendulum. It was proposed very early, but it was used to study chaos later. Because there is a very obvious and simple physical model behind it, experiments can even be done to observe the nonlinearity of this equation [3]. The form of the equation is:

Different from the previous two equations, duffing equation has a forced vibration term with time t, so it does not belong to autonomous system. It can be seen that although the system is second-order, it still has very complex nonlinearity.

As shown in the figure below, fix the amplitude frequency r and w of the excitation and change the damping d.

It can be seen that with the increase of damping d, the system changes from chaos to 2 periods and then to single period motion.

For this chaotic phenomenon, only observing the trajectory diagram can not see any law. In the next chapter, we will introduce a new observation method - Poincare section method.

Code attached:

clcclearclose all%% Lorenz attractor  h=1e-3;x0=0:h:40;[y1,~]=ODE_RK4_hyh(x0,h,[1;4;20],{'Lorenz',[10,8/3,20]});Lx=y1(1,:);Ly=y1(2,:);Lz=y1(3,:);
figure(1)subplot(1,3,1)plot(Lx,Ly);title('r=20')figure(2)subplot(1,3,1)plot3(Lx,Ly,Lz);view([51,30]);title('r=20')
[y1,~]=ODE_RK4_hyh(x0,h,[-13;-2;41],{'Lorenz',[10,8/3,28]});Lx=y1(1,:);Ly=y1(2,:);Lz=y1(3,:);
figure(1)subplot(1,3,2)plot(Lx,Ly);title('r=28')figure(2)subplot(1,3,2)plot3(Lx,Ly,Lz);view([51,30]);title('r=28')
[y1,~]=ODE_RK4_hyh(x0,h,[1;4;67],{'Lorenz',[10,8/3,99.36]});Lx=y1(1,:);Ly=y1(2,:);Lz=y1(3,:);
figure(1)subplot(1,3,3)plot(Lx,Ly);title('r=99.36')figure(2)subplot(1,3,3)plot3(Lx,Ly,Lz);view([51,30]);title('r=99.36')
%% Rossler attractor  h=2e-3;x0=0:h:180;[y1,~]=ODE_RK4_hyh(x0,h,[4.9;-5;0.07],{'Rossler',[0.1,0.1,4]});Lx=y1(1,:);Ly=y1(2,:);Lz=y1(3,:);
figure(4)subplot(2,2,1)plot3(Lx,Ly,Lz);view([51,30]);title('c=4 Cycle 1')
[y1,~]=ODE_RK4_hyh(x0,h,[9.1;-5;0.17],{'Rossler',[0.1,0.1,6]});Lx=y1(1,:);Ly=y1(2,:);Lz=y1(3,:);
figure(4)subplot(2,2,2)plot3(Lx,Ly,Lz);view([51,30]);title('c=6 Cycle 2')
[y1,~]=ODE_RK4_hyh(x0,h,[12.8;-5;0.277],{'Rossler',[0.1,0.1,8.5]});Lx=y1(1,:);Ly=y1(2,:);Lz=y1(3,:);
figure(4)subplot(2,2,3)plot3(Lx,Ly,Lz);view([51,30]);title('c=8.5 Cycle 4')
[y1,~]=ODE_RK4_hyh(x0,h,[12.8;-5;0.277],{'Rossler',[0.1,0.1,9]});Lx=y1(1,:);Ly=y1(2,:);Lz=y1(3,:);
figure(4)subplot(2,2,4)plot3(Lx,Ly,Lz);view([51,30]);title('c=9 chaos')
%% Duffing attractor  h=2e-3;x0=0:h:180;[y1,~]=ODE_RK4_hyh(x0,h,[1;0.5],{'Duffing',[1.15,1,1]});Lx=y1(1,:);Ly=y1(2,:);figure(6)subplot(1,3,1)plot(Lx,Ly);title('d=1.15')
[y1,~]=ODE_RK4_hyh(x0,h,[0.8;0.75],{'Duffing',[1.35,1,1]});Lx=y1(1,:);Ly=y1(2,:);figure(6)subplot(1,3,2)plot(Lx,Ly);title('d=1.35')
[y1,~]=ODE_RK4_hyh(x0,h,[0.7;0.73],{'Duffing',[1.5,1,1]});%[0.7;0.73]Lx=y1(1,:);Ly=y1(2,:);figure(6)subplot(1,3,3)plot(Lx,Ly);title('d=1.5')
function [F,Output]=Fdydx(x,y,Input)%Form Y'=F(x,Y)For the equation, see the relevant knowledge of numerical analysis and solving constant coefficient differential equations%The higher order is represented by a column vector, F=[dy(1);dy(2)];y(1)Is a function, y(2)Is the derivative of the function switch Input{1}    case 'Lorenz'        a=Input{2}(1);b=Input{2}(2);r=Input{2}(3);        dy(1)=a*(y(2)-y(1));        dy(2)=r*y(1)-y(2)-y(1)*y(3);        dy(3)=y(1)*y(2)-b*y(3);        F=[dy(1);dy(2);dy(3)];    case 'Rossler'        a=Input{2}(1);b=Input{2}(2);c=Input{2}(3);        dy(1)=-y(2)-y(3);        dy(2)=y(1)+a*y(2);        dy(3)=b+y(3)*(y(1)-c);        F=[dy(1);dy(2);dy(3)];    case 'Duffing'        d=Input{2}(1);r=Input{2}(2);w=Input{2}(3);        dy(1)=y(2);        dy(2)=-y(1)^3+y(1)-d*y(2)+r*cos(w*x);        F=[dy(1);dy(2)];endOutput=[];end
function [y,Output]=ODE_RK4_hyh(x,h,y0,Input)%4 rank RK method%h Algorithm with constant interval y=zeros(size(y0,1),size(x,2));y(:,1)=y0;for ii=1:length(x)-1    yn=y(:,ii);    xn=x(ii);    [K1,~]=Fdydx(xn    ,yn       ,Input);    [K2,~]=Fdydx(xn+h/2,yn+h/2*K1,Input);    [K3,~]=Fdydx(xn+h/2,yn+h/2*K2,Input);    [K4,~]=Fdydx(xn+h  ,yn+h*K3  ,Input);    y(:,ii+1)=yn+h/6*(K1+2*K2+2*K3+K4);endOutput=[];end