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

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?

[darova]

Here is your equation

Why can't you just use ode45 solver? What is the main problem? Can you explain more?

[Roderick]

Hi, darova. Thank for your comment. Yes, if you look at the part of my code that is my commented with %, I have used ode45. The problem is that I want to solve that equation, in the two regimes explained in Eq. (2) and (3), for each value of pulse_duration, amplitude, and Frequency. My problem is that I do not exactly how I can store the solution of this differential equation for each combination of the aforementioned parameters.

[Roderick]

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.

Answer 1

[darova]

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?

[Roderick]

Hey! Thank you very much! That helped my a lot. Yes, Hz it is involved in the definition of mz, which is defined in the for loop. Do you suggest defining Hz in the same way that you have defined Hzdot? What would the syntax for mz be like if that were so?

[darova]

• Do you suggest defining Hz in the same way that you have defined Hzdot?

Yes, definitely

And expression for mz:

[Roderick]

Great! Everything is alright now, I guess. Thank you very much!