Solve a system of three coupled first-order differential equations

Hello eveyone,

I'm trying to solve the system of coupled first-order differential equations that appears in the PDF that I attach to this post. To do this, I have preliminarily created the following scripts:

function dxdt=myode(t,x,alpha,h,a,lambda)
    dxdt(1)=(x(3)/(1+alpha^2))*(sin(2*x(2))/2+alpha*h);
    dxdt(2)=(1/(1+alpha^2))*(h-(alpha/2)*sin(2*x(2)));
    dxdt(3)=(1/(alpha*a))*(1/x(3)-(1+lambda*sin(2*x(2)))*x(3));
    dxdt=dxdt(:);
end

clc;
clear all;
close all force;
%% First of all, let's introduce the related parameters
alpha=0.001;
lambda=[1 10 100 1000 10000];
a=(pi^2)/6;
for i=1:length(lambda)
   field(:,i)=[alpha*lambda(i)/4,alpha*lambda(i)/2,alpha*lambda(i),3*alpha*lambda(i)/2]'; 
end
lambda_name={'1' '10' '100' '1000' '10000'};
field_name={'alpha_lambda_4' 'alpha_lambda_2' 'alpha_lambda' '3_alpha_lambda_2'};
%% Now, let's introduce the related time array
time=linspace(0,5000,10000);
%% Now, let's introduce the initial conditions
ic=[0,0,1];
%% Now, let's solve the coupled first-order differential equations
for i=1:length(lambda)
    for j=1:length(field(:,i))
        [t,sol]=ode45(@(t,x)myode(t,x,alpha,field(j,i),a,lambda(i)),time,ic);
        fnm=sprintf('C:\\Users\\890384\\Documents\\MEGA\\PhD\\Papers\\Micromagnetic Cell Size\\Numerical Calculations Thiaville Approach\\Thiaville_Numerical_lambda=%s_field=%s.txt',lambda_name{i},field_name{j});
        dlmwrite(fnm,[time' sol],'delimiter',' ')
    end
end

However, I don't know if it makes sense how I have written my system of differential equations in the function domain. There, I am not specifying, for example, that dxdt(1) is the first derivative with respect to time t of x(1), and the same about dxdt(2), x(2), dxdt(3), and x(3). To write the first two differential equations dxdt(1) and dxdt(2) I have decoupled the first two equations that appear in the PDF. The program works and yields results. But there are some cases of the ones contemplated in my script that output normal numerical results first, and then NaN thereafter.

Could you give me some feedback on this? If it wasn't right, how could you reorient it?

Answer 1

[Alan Stevens]

According to the pdf equations, instead of

dxdt(3)=(1/(alpha*a))*(1/x(3)-(1+lambda*sin(2*x(2)))*x(3));

you should have

 dxdt(3)=(1/(alpha*a))*(1/x(3)-(1+lambda*sin(x(2))^2)*x(3));

Also, I suggest you use ode15s rather than ode45.

[Richard Wood]

Hi Alan, thank you for your comments. You are right, I have corrected the dxdt(3) equation, thanks for noticing. On the other hand, if you allow me, I have some questions related to your answer. Why do you advise me to use ode15s instead of ode45 in this case? And on the other hand, it continues to surprise me that my differential equations are well expressed in the function environment without specifying that dxdt(1) is the first derivative in time of x(1) and so on, is my approach equally good?

[Alan Stevens]

I suggested ode15s because ode45 seemed very much slower for the higher values of lambda.

Also, I think you overdo the time specification. Why do you need an output of 10000 timesteps?

Other than that, your code seems to work (though I didn't test the two lines that print the results).