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

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?

Answer 1

[Ameer Hamza]

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

[Richard Wood]

Hey! Thank you very much! I had an idea on how to deal with it, but I was a little bit lost at the time to implement it. I appreciate your words!

[Ameer Hamza]

I am glad to be of help.