Solve a system of three coupled first-order differential equations

18 次查看（过去 30 天）
Roderick
Roderick 2020-11-12
评论： Alan Stevens 2020-11-13
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
Alan Stevens 2020-11-12
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 个评论
Roderick
Roderick 2020-11-13
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 2020-11-13
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).

请先登录，再进行评论。

更多回答（0 个）

类别

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

标签

Community Treasure Hunt

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

Start Hunting!

Translated by