Solve a system of three coupled first-order differential equations

35 visualizaciones (últimos 30 días)
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?

Respuesta aceptada

Alan Stevens
Alan Stevens el 12 de Nov. de 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 comentarios
Richard Wood
Richard Wood el 13 de Nov. de 2020
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
Alan Stevens el 13 de Nov. de 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).

Iniciar sesión para comentar.

Más respuestas (0)

Categorías

Más información sobre Programming en Help Center y File Exchange.

Community Treasure Hunt

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

Start Hunting!

Translated by