# Solve a system of three coupled first-order differential equations

51 views (last 30 days)
Richard Wood on 12 Nov 2020
Commented: Alan Stevens on 13 Nov 2020
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?

Alan Stevens on 12 Nov 2020
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.
##### 2 CommentsShowHide 1 older comment
Alan Stevens on 13 Nov 2020
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).