**You are now following this question**

- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.

# Phase diagram of a second-order differential equation

11 views (last 30 days)

Show older comments

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?

##### 6 Comments

Richard Wood
on 11 Mar 2020

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
on 11 Mar 2020

- 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

Richard Wood
on 11 Mar 2020

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
on 11 Mar 2020

- 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

Richard Wood
on 11 Mar 2020

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.

### Accepted Answer

darova
on 11 Mar 2020

Edited: darova
on 11 Mar 2020

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

##### 9 Comments

Richard Wood
on 12 Mar 2020

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
on 12 Mar 2020

- 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

Richard Wood
on 12 Mar 2020

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
on 12 Mar 2020

- 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

Richard Wood
on 12 Mar 2020

Richard Wood
on 12 Mar 2020

darova
on 12 Mar 2020

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

Richard Wood
on 12 Mar 2020

### More Answers (0)

### See Also

### Categories

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!**An Error Occurred**

Unable to complete the action because of changes made to the page. Reload the page to see its updated state.

Select a Web Site

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

You can also select a web site from the following list:

## How to Get Best Site Performance

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

### Americas

- América Latina (Español)
- Canada (English)
- United States (English)

### Europe

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)

### Asia Pacific

- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)