Phase diagram of a second-order differential equation

Hello everyone

I have solved a second-order differential equation, and as a result of it I have obtained the values of an angle, phi, and its first derivative on time, phidot, assuming that a time equal to zero both are zero. Now, I would like to do a phase diagram as the one that I have attached. Which is the most suitable function to plot and what I need?

Any suggestion?

[darova]

What about quiver and streamline?

See my answer here

[Roderick]

Hey! Thank you for your comment. Definitely quiver is what I think I need, as I imagined. The thing is that, as the user that asked you that question, I am not very sure of what each parameter on the syntax of quiver means (I mean, u, v, x, and y). Imagine that I have already the angle (theta, as in the figure that I posted) and its first derivative (theta dot, as in the figure that I posted) coming from the solution of a second-order differential equation. Each of them is given by an array of (70001x1 elements). Which is the role of this variables that I have already calculated on the syntax of the quiver function? Should I create a meshgrid beginning from this two arrays encapsulating the solutions of the original second-order differential equation?

What is the meaning of DX and DY in your example in the other post?

[darova]

• I am not very sure of what each parameter on the syntax of quiver means (I mean, u, v, x, and y).

x and y are positions of each vector

u and v are vectors

[Roderick]

Thank for your answer. Okay, so let's focus on the syntax that you have proposed on the aforementioned post:

[X,Y] = meshgrid(-2:0.2:2);
F = @(X,Y) 4*Y.^2 - 16;
[x,y] = ode45(F,[-2 2],2-1e-4);
DY = F(1,Y);
DX = DY*0+1;

I assume that first I have to create the meshgrid of which is the complete space on which I want to map my solution. So, let's say, I could create a meshgrid going from the minimum of the angular variable to its maximum value for the horizontal axis, and the same on the vertical axis with its first derivative. This is only plausible if the (u,v) coordinates of an arrow does not go outside of this meshgrid, so probably it is neccesary to create one more extensive. After that, you define a variable F, which depends on the previous meshgrid, and where you are going to evaluate an expression depending only on Y (maybe Y depends on X somehow). Once you do that, you solve by ode45, putting as startpoints? [-2 2], which coincides with the meshgrid limits, and then you specify 2-1e-4, which could some kind of initial condition. After that, you define DY, evaluating F on the X=1 point for all the values created on the meshgrid on Y, and then you define DX somehow (I do not understand that point at all). Up to which point do you that this approach is compatible with the variables that I have already? Could I use this script directly taking into account that I have already the values of the angular variable and its derivative already calculated? Sorry for bothering you with silly questions.

[darova]

• After that, you define a variable F, which depends on the previous meshgrid, and where you are going to evaluate an expression depending only on Y (maybe Y depends on X somehow)

Variable F in this case is just a derivative (depends on y only, but should be 2 inputs for ode45)

• After that, you define DY, evaluating F on the X=1 point for all the values created on the meshgrid on Y

The correct form would be as following (i think)

DY = F([],Y);               % doesn't matter what X=1 (equation doesn't have X variable)

• and then you define DX somehow (I do not understand that point at all).

I evaluated , but to plot quiver i need u and v (dx and dy)

So i assign dx=1(always) and dy=4y-16

• Up to which point do you that this approach is compatible with the variables that I have already? Could I use this script directly taking into account that I have already the values of the angular variable and its derivative already calculated?

please show your equation

[Roderick]

Thank you very much for your comment. I think that it could be better for our mutual understanding if I provide you with as much information of my real problem as possible. First of all, I have attached some data, Example.m. There, the first column corresponds to time (the variable with the respect to which the derivatives are performed), the second column corresponds to the angular values and the third one to its time derivative. This values has been calculated from the following equation:

The values involved are the following:

Ss=5/2;
mu0=4*pi*10^(-7); %H/m
gamma=2.21*10^5; %m/(A*s)
kB=1.38064852*10^(-23); %J/K
a0=3.328*10^(-10); %m
c=8.539*10^(-10); %m
V=c*(a0^2); %m^3
Hx=0;
Hy=40.*(40*79.5774715459)); % A/m
muB=9.27400994*10^(-24); %(T^2*m^3)/J
mu=4*muB; %(T^2*m^3)/J
Ms=mu/V; %T^2/J
K4parallel=1.8548*(10^(-25))/V; %J/m^3
alpha=0.001;
Omegae=(2*gamma*abs(4*(-396*kB/V)-532*kB/V))/(mu0*Ms); %s^-1
Omega4parallel=(2*gamma*K4parallel)/(mu0*Ms); %s^-1
OmegaR=sqrt(2*Omegae*Omega4parallel);

Maybe there are some parameters that are not needed, sorry if that it is the case. So at the end, this is what I have to face the quiver plot.

Answer 1

[darova]

It's your equation

Assume you have a curve

Then to create a quiver or streamline you need ( ϕ, , u, v )

I looked here and see that (ϕ, ) in [ ]

What is the connection between u, v and ?

So substituting (ϕ, ) in we have angle

I assume u=1, then v=a

I used parameters you gave

F = @(phi,dphi) -OmegaR^2/4*sin(4*phi) -2*Omegae*gamma*Hy*cos(phi) -2*Omegae*alpha*dphi;
xx = -pi:0.3:pi;
[p,dp] = meshgrid(xx);      % grid for phi and dphi
v = F(p,dp)./dp;
u = v*0 + 1;
quiver(p,dp,u,v,'b')
hold on
streamline(p,dp,u,v,xx,xx,'r')
hold off

i got this

streamline somehow didn't work

[Roderick]

Thank you very much for your complete answer! However, I have made a mistake in the data I gave you earlier. Specifically in Hy. To resolve the second-order differential equation I have created an array of times to specify which is the profile of Hy on time, and then use it to solve the equation:

time=[0:0.001:70].*(10^(-12)); % s
critical_index_time=find(time==3E-12);
time_first_interval=time(1:critical_index_time); % s
time_second_interval=time((critical_index_time+1):end); % s
Hy1=(ones(1,length(time_first_interval)).*(40*79.5774715459)); % A/m
Hy2=zeros(1,length(time_second_interval));
Hy=[Hy1 Hy2]';

Do you think that it is straightforward to implement this change on your previous approach?

Moreover, on which point of your code do you have used my data of phi and dphi? Sorry for the silly question.

[darova]

• To resolve the second-order differential equation I have created an array of times to specify which is the profile of Hy on time, and then use it to solve the equation:

Hm. That changes everything. I think ode45 should be used for each curve/solution

• Moreover, on which point of your code do you have used my data of phi and dphi?

Don't understand the question. Are you asking where i used phi and dphi?

xx = -pi:0.3:pi;
[p,dp] = meshgrid(xx);      % grid for phi and dphi

[Roderick]

If helps, my actual array looks like this:

clc
clear all
close all force
outputdir='Figures';
%% First, let's write which are the initial values of our problem
% For t=0, we have that phi=x(1)=0, and that \dot{phi}=x(2)=0
xinitial=0;
dxinitial=0;
initial_conditions1=[xinitial dxinitial];
% Also, I need to define a parameter, Omegae
Neel_Temperature=1575; %K
beta=0.333;
Ss=5/2;
mu0=4*pi*10^(-7); %H/m
gamma=2.21*10^5; %m/(A*s)
kB=1.38064852*10^(-23); %J/K
a0=3.328*10^(-10); %m
c=8.539*10^(-10); %m
V=c*(a0^2); %m^3
Hx=0;
muB=9.27400994*10^(-24); %(T^2*m^3)/J
mu=4*muB; %(T^2*m^3)/J
Ms=mu/V; %T^2/J
K4parallel=1.8548*(10^(-25))/V; %J/m^3
alpha_T0=0.001;
Omegae_T0=(2*gamma*abs(4*(-396*kB/V)-532*kB/V))/(mu0*Ms); %s^-1
Omega4parallel_T0=(2*gamma*K4parallel)/(mu0*Ms); %s^-1
OmegaR_T0=sqrt(2*Omegae_T0*Omega4parallel_T0);
phipotential=[0:0.01:pi/2];
thetapotential=pi/2;
UwithoutDW=-(Ss^2)*((K4parallel*(Ss^2)/2)*(sin(thetapotential).*sin(phipotential)).^4+...
    (K4parallel*(Ss^2)/2)*(sin(thetapotential).*cos(phipotential)).^4)*V;
%% Now, let's calculate the parameters that depend on temperature
Temperature=[0:200:1000]; %K
magnetization_temperature=(1-(Temperature./Neel_Temperature)).^(beta);
alpha=alpha_T0.*(1-(Temperature./(3*Neel_Temperature)));
Omegae=Omegae_T0.*magnetization_temperature; % s^(-1)
OmegaR=OmegaR_T0.*(magnetization_temperature.^5); % s^(-1)
% Temperature potential
UTemperature_200_K=ones(length(phipotential),1).*(kB*Temperature(2))+0.*phipotential+UwithoutDW(1);
UTemperature_400_K=ones(length(phipotential),1).*(kB*Temperature(3))+0.*phipotential+UwithoutDW(1);
UTemperature_600_K=ones(length(phipotential),1).*(kB*Temperature(4))+0.*phipotential+UwithoutDW(1);
UTemperature_800_K=ones(length(phipotential),1).*(kB*Temperature(5))+0.*phipotential+UwithoutDW(1);
UTemperature_1000_K=ones(length(phipotential),1).*(kB*Temperature(6))+0.*phipotential+UwithoutDW(1);
%% Now, let's create the array of times on which we are interested
time=[0:0.001:70].*(10^(-12)); % s
critical_index_time=find(time==3E-12);
time_first_interval=time(1:critical_index_time); % s
time_second_interval=time((critical_index_time+1):end); % s
Hy1=(ones(1,length(time_first_interval)).*(40*79.5774715459)); % A/m
Hy2=zeros(1,length(time_second_interval));
Hy=[Hy1 Hy2]';
%% Let's solve the second-order differential equation, where the X variable will have to columns
% Probably, the first column corresponds to phi=x(1) and the second one to
% \dot{phi}=x(2)
for i=1:length(Temperature)
    [T1{i},X1{i}]=ode45(@(t,x)fx(t,x,Hx,Hy1(i),OmegaR(i),Omegae(i),gamma,alpha(i)),time_first_interval,initial_conditions1);
end
Xmat1=cell2mat(X1);
Tmat1=cell2mat(T1);
for i=1:length(Temperature)
    [T2{i},X2{i}]=ode45(@(t,x)fx(t,x,Hx,Hy2(i),OmegaR(i),Omegae(i),gamma,alpha(i)),time_second_interval,[Xmat1(end,2*i-1) Xmat1(end,2*i)]);
end
Xmat2=cell2mat(X2);
Tmat2=cell2mat(T2);
Xmat=vertcat(Xmat1,Xmat2);
Tmat=vertcat(Tmat1,Tmat2);
lx=cos(Xmat(:,1:2:end));
ly=sin(Xmat(:,1:2:end));
mz=-(1./(2*Omegae)).*Xmat(:,2:2:end);
phi=Xmat(:,1:2:end);
phidot=Xmat(:,2:2:end);
% fid_example=fopen('Example.m','wt');
% fprintf(fid_example, '%d %d %d \n',[Tmat(:,1) phi(:,1) phidot(:,1)]');
% fclose(fid_example);
phiddot=diffxy(Tmat,Xmat(:,2:2:end));

My idea was to do that phase diagram for each row of phi and phidot, being each row dictated by a value of the Temperature vector. So at the end, I want to have six phase diagrams.

• Don't understand the question. Are you asking where i used phi and dphi?

Yeah, exactly. Because when you define your F, you put as variables phi and dphi, but in principle you are not substituting there the values of the solutions, because after that you only evaluate F for the meshgrid that you have created by using p and dp.

[darova]

• you are not substituting there the values of the solutions

Yes. I just assumed that phi and dphi somewhere in [-pi .. pi]

Do you already have solutions? Im starting to get lost

[Roderick]

Yeah, right. I have solved the second-order differential equation using ode45, so I have the value of phi and dphi for my parameters and initial conditions.

[darova]

So all you have to do is to plot arrows?

[Roderick]

I think so, yes. Maybe the arrows and the contour lines. If you see my whole code, I want to build a for loop in i=1:length(Temperature) (which is the same as the number of columns of phi and phidot), and replicate the kind of plot that I put when I posted my answer. Sorry if I was not clear enough.

[darova]

Try this

x = 0:0.1:10;
y = sin(x);
u = diff(x);
v = diff(y);
ii = 1:10:length(u);                % how many arrows
plot(x,y)
hold on
quiver(x(ii),y(ii),u(ii),v(ii))
hold off

[Roderick]

Thank you again for your comment. But I think we are not understanding each other, surely because of me. In your last code, where is the role of my solution of phy and phidot (or dphi, as you called it at the beginning)? I suppose that you create your y function to play the role of phidot, and so x to be my phi variable. It is not neccesary for you to answer me again, probably I am wasting your time. If you run my code, eliminating the last line and the line of the outputdir you will see exactly what I have, how many elements on each array of phi and phidot, and the number of columns being the same as the number of elements of the Temperature variable.