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

3 次查看（过去 30 天）
Roderick
Roderick 2020-5-19
评论： Roderick 2020-5-20
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 个评论
显示 1更早的评论
Roderick
Roderick 2020-5-19
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
Roderick 2020-5-20
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.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

darova
darova 2020-5-20
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 个评论
显示 1更早的评论
darova
darova 2020-5-20
  • Do you suggest defining Hz in the same way that you have defined Hzdot?
Yes, definitely
And expression for mz:

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Mathematics 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by