Solve a system of three coupled first-order differential equations
16 views (last 30 days)
Show older comments
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?
0 Comments
Accepted Answer
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 Comments
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).
More Answers (0)
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!