ode45 inside a nested loop

7 visualizaciones (últimos 30 días)
Richard Wood
Richard Wood el 28 de Feb. de 2020
Comentada: Richard Wood el 29 de Feb. de 2020
Hello everyone,
My situation is the following. I want to solve a second-order differential equation, where the parameters involved depends on external parameters (in my case, on variables Temperature and Hy1_value). This parameters are Omegae, OmegaR, alpha. My function file looks as follows
function dx=fx(t,x,Hx,Hy,OmegaR,Omegae,gamma,alpha);
% Here, x represents the in-plane angle, phi, and t represents the
% time variable
% Global variables: Hx, Hy, OmegaR, Omegae, gamma, alpha
% Here, x(1)=phi, x(2)=\dot{phi}
dx=zeros(2,1);
dx(1)=x(2);
dx(2)=-(OmegaR.^2/4).*sin(4*x(1))+2*Omegae.*gamma.*(Hy.*cos(x(1))-Hx.*sin(x(1)))-...
2*Omegae.*alpha.*x(2);
end
Now, I want to solve my second-order differential equation which is inside a nested for loop. The thing is that I am not sure of how to index the variables involved in the second-order differential equation solving process by using ode45. My code looks as follows:
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;
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);
pulse_duration=1E-12; % s
decay_time=100E-12; % s
time=[0:0.001:70].*(10^(-12)); % s
critical_index_time=find(time==1E-12);
time_first_interval=time(1:critical_index_time); % s
time_second_interval=time((critical_index_time+1):end); % s
Temperature_max=[200:200:1200]; % K
Hy1_values=linspace(100,1000,10).*79.5774715459; % A/m
Temperature1=zeros(length(time_first_interval),length(Temperature_max));
Temperature2=zeros(length(time_second_interval),length(Temperature_max));
Temperature=zeros(length(time),length(Temperature_max));
Hy1=zeros(length(time_first_interval),length(Hy1_values),length(Temperature_max));
magnetization_temperature=zeros(length(time),length(Temperature_max));
magnetization_temperature1=zeros(length(time_first_interval),length(Temperature_max));
magnetization_temperature2=zeros(length(time_second_interval),length(Temperature_max));
alpha=zeros(length(time),length(Temperature_max));
Omegae=zeros(length(time),length(Temperature_max)); % s^(-1)
OmegaR=zeros(length(time),length(Temperature_max)); % s^(-1)
alpha1=zeros(length(time_first_interval),length(Temperature_max));
Omegae1=zeros(length(time_first_interval),length(Temperature_max)); % s^(-1)
OmegaR1=zeros(length(time_first_interval),length(Temperature_max)); % s^(-1)
alpha2=zeros(length(time_second_interval),length(Temperature_max));
Omegae2=zeros(length(time_second_interval),length(Temperature_max)); % s^(-1)
OmegaR2=zeros(length(time_second_interval),length(Temperature_max)); % s^(-1)
Hy1=zeros(length(time_first_interval),length(Hy1_values),length(Temperature_max)); % A/m
Hy2=zeros(length(time_second_interval),length(Hy1_values),length(Temperature_max));
Xmat1=zeros(length(time_first_interval),2*length(Hy1_values),length(Temperature_max));
%% Now, let's calculate the parameters that depend on temperature
for i=1:length(Temperature_max)
for j=1:length(Hy1_values)
Temperature1(:,i)=(Temperature_max(i)./sqrt(pulse_duration)).*sqrt(time_first_interval); % K
Temperature2(:,i)=Temperature1(end,i).*(exp(-time_second_interval./pulse_duration)); % K
Temperature(:,i)=vertcat(Temperature1(:,i),Temperature2(:,i));% K
magnetization_temperature1(:,i)=(1-(Temperature1(:,i)./Neel_Temperature)).^(beta);
magnetization_temperature2(:,i)=(1-(Temperature2(:,i)./Neel_Temperature)).^(beta);
magnetization_temperature(:,i)=vertcat(magnetization_temperature1(:,i),...
magnetization_temperature2(:,i));
alpha1(:,i)=alpha_T0.*(1-(Temperature1(:,i)./(3*Neel_Temperature)));
alpha2(:,i)=alpha_T0.*(1-(Temperature2(:,i)./(3*Neel_Temperature)));
alpha(:,i)=vertcat(alpha1(:,i),alpha2(:,i));
Omegae1(:,i)=Omegae_T0.*magnetization_temperature1(:,i); % s^(-1)
Omegae2(:,i)=Omegae_T0.*magnetization_temperature2(:,i); % s^(-1)
Omegae(:,i)=vertcat(Omegae1(:,i),Omegae2(:,i)); % s^(-1)
OmegaR1(:,i)=OmegaR_T0.*(magnetization_temperature1(:,i).^5); % s^(-1)
OmegaR2(:,i)=OmegaR_T0.*(magnetization_temperature2(:,i).^5); % s^(-1)
OmegaR(:,i)=vertcat(OmegaR1(:,i),OmegaR2(:,i));
Hy1(:,j,i)=ones(length(time_first_interval),1).*Hy1_values(j); % A/m
Hy2(:,j,i)=zeros(length(time_second_interval),1); % A/m
Hy(:,j,i)=vertcat(Hy1(:,j,i),Hy2(:,j,i));
%% 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)
% I DON'T KNOW HOW TO DEAL WITH THE INDEXING FROM HERE ON!
[T1{i,j},X1{i,j}]=ode45(@(t,x)fx(t,x,Hx,Hy1(:,j,i),OmegaR1(:,i),Omegae1(:,i),gamma,alpha1(:,i)),time_first_interval,initial_conditions1);
% Xmat1(:,2*j-1,i)=cell2mat(X1);
% Tmat1=cell2mat(T1);
% [T2,X2]=ode45(@(t,x)fx(t,x,Hx,Hy2,OmegaR,Omegae,gamma,alpha),time_second_interval,[Xmat1(end,1) Xmat1(end,2)]);
% Xmat2=cell2mat(X2);
% Tmat2=cell2mat(T2);
% Xmat=vertcat(Xmat1,Xmat2);
% Tmat=vertcat(Tmat1,Tmat2);
% lx=cos(Xmat(:,1));
% ly=sin(Xmat(:,1));
% mz=-(1./(2*Omegae)).*Xmat(:,2);
end
end
I had to separate the solving process in two ode45's because of the abrupt discontinuity imposed in the Hy variable in the intersection of Hy1 and Hy2. The thing is that all variables are well defined. The problem begins when I put the warning "% I DON'T KNOW HOW TO DEAL WITH THE INDEXING FROM HERE ON!". My objetive is that, at the end, my arrays lx, ly and mz have the following dimensions:
lx=zeros(length(time),2*length(Hy1_values),length(Temperature_max));
ly=zeros(length(time),2*length(Hy1_values),length(Temperature_max));
mz=zeros(length(time),2*length(Hy1_values),length(Temperature_max));
So that I have all the information after the for loops well-organized and saved.
Any suggestions?
The error that MatLab gives to me is
Unable to perform assignment because the left and right sides have a different number of elements.
Error in fx (line 8)
dx(2)=-(OmegaR.^2/4).*sin(4*x(1))+2*Omegae.*gamma.*(Hy.*cos(x(1))-Hx.*sin(x(1)))-...
Error in
Second_Order_Differential_Equation_Laser_Like>@(t,x)fx(t,x,Hx,Hy1(:,j,i),OmegaR1(:,i),Omegae1(:,i),gamma,alpha1(:,i))
(line 80)
[T1{i,j},X1{i,j}]=ode45(@(t,x)fx(t,x,Hx,Hy1(:,j,i),OmegaR1(:,i),Omegae1(:,i),...
Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in Second_Order_Differential_Equation_Laser_Like (line 80)
[T1{i,j},X1{i,j}]=ode45(@(t,x)fx(t,x,Hx,Hy1(:,j,i),OmegaR1(:,i),Omegae1(:,i),...
The thing is that OmegaR, Omegae and alpha depends on time. I have seen that when you are dealing time-dependence of the parameters involved on your differential equation, when using ode45, you must put (t) after then to reflect its dependence on time, but I think that I should also take care of the indexing of that multidimensional arrays!
  2 comentarios
darova
darova el 28 de Feb. de 2020
Here is the problem
I suppose there is some dependence between Hx and t (or x)?
Can you show original equations/formulas,
I also suggest you to give shorter names to your variables. Two weeks will pass and you will not be able to read your code
t1int = time(1:critical_index_time); % s
t2int = time((critical_index_time+1):end); % s
Temp_max=[200:200:1200]; % K
[Temp1, mag_temp1, alpha1, Omegae1, OmegaR1] = deal( zeros(length(t1int),length(Temp_max)) );
[Temp2, mag_temp2, alpha2, Omegae2, OmegaR2] = deal( zeros(length(t2int),length(Temp_max)) );
[Temp, mag_temp, alpha, Omegae, OmegaR] = deal( zeros(length(time),length(Temp_max)) );
Richard Wood
Richard Wood el 28 de Feb. de 2020
Editada: Richard Wood el 28 de Feb. de 2020
Hey! Thank for your answer. I have simplified a little bit one of the scripts, as you have proposed. I hope that that helps to read it easily. Again, my function script:
function dx=fx(t,x,Hx,Hy,OmegaR,Omegae,gamma,alpha);
% Here, x represents the in-plane angle, phi, and t represents the
% time variable
% Global variables: Hx, Hy, OmegaR, Omegae, gamma, alpha
% Here, x(1)=phi, x(2)=\dot{phi}
dx=zeros(2,1);
dx(1)=x(2);
dx(2)=-(OmegaR.^2/4).*sin(4*x(1))+2*Omegae.*gamma.*(Hy.*cos(x(1))-Hx.*sin(x(1)))-...
2*Omegae.*alpha.*x(2);
end
The other script:
%% 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 certain parameters of my material
NeelT=1575; %K
beta=0.333;
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);
pulse_duration=1E-12; % s
decay_time=100E-12; % s
%% Now,I can define the time array and also the parameters that will depend
% on Temperature (Temp, T1int, T2int). That it is, OmegaR(1,2),
% Omegae(1,2), alpha(1,2). The thing is that temperature will depend on
% time, so these coefficients that depend on temperature will depend on
% time
time=[0:0.001:70].*(10^(-12)); % s
critical_index_time=find(time==1E-12);
t1int=time(1:critical_index_time)'; % s
t2int=time((critical_index_time+1):end)'; % s
Tmax=[200:200:1200]; % K
Hy1_values=linspace(100,1000,10).*79.5774715459; % A/m
Hy1=zeros(length(t1int),length(Hy1_values),length(Tmax)); % A/m
Hy2=zeros(length(t2int),length(Hy1_values),length(Tmax));
Hy=zeros(length(time),length(Hy1_values),length(Tmax));
[magT1, Omegae1, OmegaR1, alpha1, T1int]=deal(zeros(length(t1int),length(Tmax)));
[magT2, Omegae2, OmegaR2, alpha2, T2int]=deal(zeros(length(t2int),length(Tmax)));
[magT, Omegae, OmegaR, alpha, Temp]=deal(zeros(length(time),length(Tmax)));
%% Now, let's calculate the parameters that depend on temperature (indirectly on time)
for i=1:length(Tmax)
for j=1:length(Hy1_values)
T1int(:,i)=(Tmax(i)./sqrt(pulse_duration)).*sqrt(t1int);
T2int(:,i)=T1int(end,i).*(exp(-t2int./pulse_duration));
Temp(:,i)=vertcat(T1int(:,i),T2int(:,i));
magT1(:,i)=(1-(T1int(:,i)./NeelT)).^(beta);
magT2(:,i)=(1-(T2int(:,i)./NeelT)).^(beta);
magT(:,i)=vertcat(magT1(:,i),magT2(:,i));
alpha1(:,i)=alpha_T0.*(1-(T1int(:,i)./(3*NeelT)));
alpha2(:,i)=alpha_T0.*(1-(T2int(:,i)./(3*NeelT)));
alpha(:,i)=vertcat(alpha1(:,i),alpha2(:,i));
Omegae1(:,i)=Omegae_T0.*magT1(:,i);
Omegae2(:,i)=Omegae_T0.*magT2(:,i);
Omegae(:,i)=vertcat(Omegae1(:,i),Omegae2(:,i)); % s^(-1)
OmegaR1(:,i)=OmegaR_T0.*(magT1(:,i).^5); % s^(-1)
OmegaR2(:,i)=OmegaR_T0.*(magT2(:,i).^5); % s^(-1)
OmegaR(:,i)=vertcat(OmegaR1(:,i),OmegaR2(:,i));
Hy1(:,j,i)=ones(length(t1int),1).*Hy1_values(j); % A/m
Hy2(:,j,i)=zeros(length(t2int),1); % A/m
Hy(:,j,i)=cat(2,[Hy1(:,j,i);Hy2(:,j,i)]);
%% Let's solve the second-order differential equation, where the X variable will have to columns
[T1{j},X1{j}]=ode45(@(t,x)fx(t,x,Hx,Hy1,OmegaR1(:,i),Omegae1(:,i),...
gamma,alpha1(:,i)),t1int,initial_conditions1);
% % Xmat1(:,2*j-1,i)=cell2mat(X1);
% % Tmat1=cell2mat(T1);
% % [T2,X2]=ode45(@(t,x)fx(t,x,Hx,Hy2,OmegaR2,Omegae2,gamma,alpha2),t2int,[Xmat1(end,1) Xmat1(end,2)]);
% % Xmat2=cell2mat(X2);
% % Tmat2=cell2mat(T2);
% % Xmat=vertcat(Xmat1,Xmat2);
% % Tmat=vertcat(Tmat1,Tmat2);
% % lx=cos(Xmat(:,1));
% % ly=sin(Xmat(:,1));
% % mz=-(1./(2*Omegae)).*Xmat(:,2);
end
end
I am sorry, but I do not see the point on your previous kind comment. Which is exactly the problem with the indexing? In my problem Hx=0, always. I suppose that your previous questions refers to Hy. The point is the following. I want to solve that second-order differential equation for different temperature profiles (external for loop) for a variety of magnetic fields (internal for loop, Hy variable, which values are enconded in Hy1_values, because in the frontier between t1int and t2int Hy goes from a constant value to be zero abruptly, that is the reason of using two ode45 solvers between these two problematic regions). Hy depends on time, of course. I am not sure that I need to define such a difficult three-dimensional array for Hy. I only want to play with different constant values for Hy on the t1int for different temperature profiles.
At the end what I want to do is to plot lx, ly and mz. One plot for each value of Tmax, being on each of them all the possible cases for Hy1.
I attach a file with the formulas in LaTeX (with not a great format, but it is better than look into MatLab syntaxis). I hope this helps a little bit!

Iniciar sesión para comentar.

Respuesta aceptada

darova
darova el 29 de Feb. de 2020
I want to show you a simple example of how ode45 works. I wrote simple solver like ode45
function main
clc
F = @(t,x) [x(2); x(1)-2*sin(x(1))];
ts = [0 5];
x0 = [0.5 1];
[t1,x] = ode45(F,ts,x0); % standard ode45 solver
[t2, res] = new_ode45_solver(F,ts,x0); % solve equations with new solver
h1 = plot(t1,x,'r');
hold on
h2 = plot(t2,res,'b');
hold off
legend([h1(1) h2(2)],'ode45 solver','my new ode45 solver')
end
function [t, res] = new_ode45_solver(F,ts,x0)
n = 100; % number of steps
dt = (ts(2)-ts(1))/n; % time step
res = zeros(n,2);
t = linspace(ts(1),ts(2),n)';
x = x0;
for i = 1:n
res(i,:) = x;
x = x + dt*F(t(i),x)'; % standard Euler scheme/method
end
end
Look at this line
x = x + dt*F(t(i),x)'; % standard Euler scheme/method
Look here
dx(2)=-(OmegaR.^2/4).*sin(4*x(1))+2*Omegae.*gamma.*(Hy.*cos(x(1))-Hx.*sin(x(1)))-...
2*Omegae.*alpha.*x(2);
% And look this
x(2) = [1 2 3] % wring! - doesn't work
ode45 works as loop. So if you pass some parameters (Hx,Hy) they can't be arrays
  4 comentarios
darova
darova el 29 de Feb. de 2020
  • Is the discrepancy between ode45 and your approach somewhat based on the discrepancy between your grids?
It doesn't matter if you use time interval [0 time] or linspace(0,time). ode45 solver choose it's own timestep and then create time span (the same span if you specified). Look HELP about ode45
  • The thing is Hy0. What I wanted was to made Hy=cte up to time=1E-12, and from that point on time, made it 0 for the rest of solving process. I think that that it is not contemplated in your approach, right? Just to know.
Sorry. FIrst uploaded file without this. Just add into the function
Hy = (t<1E-12)*Hy;
  • What does the introduction of this variable imply?
I just wanted to color each curve in different color
k = 0;
cmap = jet( length(Hy0)*length(Tmax0) );
But later decided to color only by Hy0 changes. So you can remove it
Richard Wood
Richard Wood el 29 de Feb. de 2020
Thank you for your comments, Darova. I have accepted your answer. If you have any useful comment about the plots, let me know, if you want. I just have some doubt about how LX, LY, MZ values are stored in each loop to undertake my plots. That it is, if I should somehow implement the plots inside the for loop for Hy or if all the values for lx, ly, mz are storaged in that variable after that loop and I can plot them outside for loop of Hy but inside of the one of temperature, of course.

Iniciar sesión para comentar.

Más respuestas (0)

Categorías

Más información sobre Loops and Conditional Statements en Help Center y File Exchange.

Community Treasure Hunt

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

Start Hunting!

Translated by