Resolution of a second-order differential equation for different regions controlled by an external parameter

1 view (last 30 days)
Richard Wood on 6 May 2020
Commented: Ameer Hamza on 7 May 2020
Hello everyone
I am trying to perform a time integration to a second-order differential equation numerically. The objective is to integrate the same equation into three different regimes. All this is explained in the document that I attach to my question. I write here the parameters involved:
mu0=4*pi*10^(-7);
gamma=2.21*10^5;
kB=1.38064852*10^(-23);
a0=3.328*10^(-10);
alpha_001=0.001;
muB=9.27400994*10^(-24);
mu=4*muB;
c=8.539*10^(-10);
V=c*(a0^2);
Ms=mu/V;
Delta0=1.9777269E-08;
vmax=43.3*10^(3);
Ss=5/2;
jinter=532*kB;
Jinter=jinter*Ss^2/V;
vmax=43.3*10^(3);
HSO=60*10^4/12.54;
C1=(2*alpha*gamma*Jinter)/(mu0*Ms);
C2=(2*(gamma^2)*Jinter*Delta0)/(mu0*Ms);
Hcrit=2.72*10^4/12.54;
t=[0:0.001:140].*(10^(-12));
ramping=[30 40 50 60].*(10^(-12));
Any idea?
0 CommentsShowHide -1 older comments

Sign in to comment.

Accepted Answer

Ameer Hamza on 7 May 2020
One of the few clearly stated questions on this website, so I will answer it :D. You first need to express your 2nd order ODE into a system of 2 first-order ODEs. See an example here: https://www.mathworks.com/help/matlab/ref/ode45.html#bu3uj8b. The following code uses a for-loop to use each value of ramping and solve the ODE. It outputs the results in a cell array. The two columns of the solution represent the values of x and . Try the following code
ramping=[30 40 50 60].*(10^(-12));
t=(0:0.1:140).*(10^(-12));
ic = [0; 0];
sol = cell(numel(ramping), 1); % each cell contains the solution for a particular value of ramping
for i=1:numel(ramping)
[~, sol{i}] = ode45(@(t,x) odeFun(t,x,ramping(i)), t, ic);
end
figure;
tiledlayout('flow');
for i=1:numel(sol)
nexttile
plot(t, sol{i});
legend({'$x$', '$\dot{x}$'}, ...
'Interpreter', 'latex', ...
'FontSize', 12)
title(sprintf('Ramping=%g', ramping(i)));
end
function dxdt = odeFun(t, x, ramping_time)
mu0=4*pi*10^(-7);
alpha=0.001;
gamma=2.21*10^5;
kB=1.38064852*10^(-23);
a0=3.328*10^(-10);
muB=9.27400994*10^(-24);
mu=4*muB;
c=8.539*10^(-10);
V=c*(a0^2);
Ms=mu/V;
Delta0=1.9777269E-08;
Ss=5/2;
jinter=532*kB;
Jinter=jinter*Ss^2/V;
vmax=43.3*10^(3);
HSO=60*10^4/12.54;
C1=(2*alpha*gamma*Jinter)/(mu0*Ms);
C2=(2*(gamma^2)*Jinter*Delta0)/(mu0*Ms);
Hcrit=2.72*10^4/12.54;
if t < ramping_time
y = max(HSO/ramping_time*t - Hcrit, 0);
elseif t < 100e-12
y = HSO;
else
y = 0;
end
dxdt(1) = x(2);
dxdt(2) = -C1*x(2) + C2*y*sqrt(1-(x(2)/vmax)^2);
dxdt = dxdt(:);
end
2 CommentsShowHide 1 older comment
Ameer Hamza on 7 May 2020
I am glad to be of help.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!