Resolution of a second-order differential equation for multiple variable parameters

3 views (last 30 days)
Richard Wood
Richard Wood on 19 May 2020
Commented: Richard Wood on 20 May 2020
Hello everyone,
I am trying to solve a second-order differential equation, given by Eq. (1) in the document attached in this post. See the details there of the system and variables used. My codes right now looks like:
function main
%% 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_conditions=[xinitial dxinitial];
% Let's define the involved parameters
mu0=4*pi*10^(-7);
gamma=2.21*10^5;
kB=1.38064852*10^(-23);
a0=3.328*10^(-10);
c=8.539*10^(-10);
V=c*(a0^2);
muB=9.27400994*10^(-24);
mu=4*muB;
Ms=mu/V;
K4parallel=1.8548*(10^(-25))/V;
Jeff=abs(-4*396-532)*kB/V;
alpha=0.001;
Omegae=(2*gamma*Jeff)/(mu0*Ms);
Omega4parallel=(2*gamma*K4parallel)/(mu0*Ms);
OmegaR=sqrt(2*Omegae*Omega4parallel);
%% Now, let's create the array of times on which we are interested
time=[0:0.00001:100].*(10^(-12));
pulse_duration=[1E-15 10E-15 0.1E-12 1E-12 10E-12];
amplitude_field=[0.1 1 10 20 40 60 80 100].*(795.7747154822216);
Frequency=[100E9 1E12 100E12 500E12 1E15];
pulse_units=["fs" "fs" "fs" "ps" "ps"];
pulse_number=[1 10 100 1 10];
amplitude_number=[0.1 1 10 20 40 60 80 100];
Frequency_units=["GHz" "THz" "THz" "THz" "PHz"];
Frequency_number=[100 1 100 500 1];
for i=1:length(pulse_duration)
for j=1:length(amplitude_field)
for k=1:length(Frequency)
Hz=@(t,pulse_duration,amplitude_field,Frequency)(t<=pulse_duration(i))...
.*amplitude_field(j).*cos(Frequency(k).*t)+(t>pulse_duration(i)).*0;
Hzdot=@(t,pulse_duration,amplitude_field,Frequency)(t<=pulse_duration(i))...
.*(-Frequency(k).*amplitude_field(j).*sin(Frequency(k).*t))+...
(t>pulse_duration(i)).*0;
[t,x]=ode45(@(t,x)myode(t,x,Hzdot),time,initial_conditions);
lx=cos(x(:,1));
ly=sin(x(:,1));
TIME=t;
mz=-(1./(2.*Omegae)).*(x(:,2)-gamma.*Hz(i));
fnm=sprintf('Data_Pulse_%g_%s_Amplitude_%g_mT_Frequency_%g_%s.dat',pulse_number(i),pulse_units(i),amplitude_number(j),Frequency_number(k),Frequency_units(k));
mat=[TIME(:) lx(:) ly(:) mz(:)];
dlmwrite(fnm,mat,'delimiter',' ')
end
end
end
function dy=myode(t,x,Hzdot)
dy(1,1)=x(2);
dy(2,1)=-(OmegaR^2/4)*sin(4*x(1))-2*Omegae*alpha*x(2)+...
gamma*Hzdot;
end
end
In principle, I get the following error:
Undefined operator '*' for input arguments of type 'function_handle'.
Error in main/myode (line 55)
gamma*Hzdot;
Error in main>@(t,x)myode(t,x,Hzdot) (line 41)
[t,x]=ode45(@(t,x)myode(t,x,Hzdot),time,initial_conditions);
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 main (line 41)
[t,x]=ode45(@(t,x)myode(t,x,Hzdot),time,initial_conditions);
Any idea on what it is happening? Moreover, I have defined some array made of characters (namely, pulse_units and Frequency_units). With this I was willing to introduce this characters in the name of the stored file on each for loop, and I have used for that %s. Is that the correct? Will the quotation marks, "", appear in the final file name?
  3 Comments
Richard Wood
Richard Wood on 20 May 2020
Dear Darova, I have modified the original post and the code involved, and I have highlighted the main problems right now. Please, take a look at it if you want. Sorry if the previous questions was unclear. I hope that now that my code looks better will be easy to understand it.

Sign in to comment.

Accepted Answer

darova
darova on 20 May 2020
Because your Hzdot function depends on t variable
dy(2,1) = -(OmegaR^2/4)*sin(4*x(1))-2*Omegae*alpha*x(2)+gamma*Hzdot(t);
Change your Hzdot function as:
Hzdot=@(t) (t<=pulse_duration(i))*(-Frequency(k).*amplitude_field(j).*sin(Frequency(k).*t));
[t,x]=ode45(@(t,x)myode(t,x,Hzdot),time,initial_conditions);
I don't you are using Hz function. Do you need it?
  3 Comments

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!

Translated by