Phase diagram of a second-order differential equation

Question

Roderick 2020-3-10

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/510119-phase-diagram-of-a-second-order-differential-equation

评论： Roderick 2020-3-12

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 个评论
显示 4更早的评论隐藏 4更早的评论

Roderick 2020-3-11

在 MATLAB Online 中打开

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 2020-3-11

在 MATLAB Online 中打开

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 2020-3-11

编辑：Roderick 2020-3-11

在 MATLAB Online 中打开

Example.m

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 2020-3-11

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/510119-phase-diagram-of-a-second-order-differential-equation#answer_419706

编辑：darova 2020-3-11

在 MATLAB Online 中打开

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 个评论
显示 7更早的评论隐藏 7更早的评论

Roderick 2020-3-12

在 MATLAB Online 中打开

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 2020-3-12

在 MATLAB Online 中打开

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 2020-3-12

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.

请先登录，再进行评论。

Phase diagram of a second-order differential equation

6 个评论
显示 4更早的评论隐藏 4更早的评论

采纳的回答

9 个评论
显示 7更早的评论隐藏 7更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

Phase diagram of a second-order differential equation

6 个评论 显示 4更早的评论隐藏 4更早的评论

采纳的回答

9 个评论 显示 7更早的评论隐藏 7更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

6 个评论
显示 4更早的评论隐藏 4更早的评论

9 个评论
显示 7更早的评论隐藏 7更早的评论