Solve a second-order differential equation with an array parameter that does not depend explicitly on the independent variable

Hello everyone, I am trying to solve a second-order differential equation, which includes a parameter whose value for each step of the independent variable was extracted previously through simulations. In my code, that variable is called Delta, which have the same length that the time variable, which is the aforementioned independent variable, which is given also by the simulation step. The function where I have defined my second-order differential equation is the following:

function dxdt = odeFun(t, x, Delta)
    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;
    Ss=5/2;
    jinter=532*kB;
    Jinter=jinter*Ss^2/V;
    C1=(2*alpha*gamma*Jinter)/(mu0*Ms);
    C2=(2*(gamma^2)*Jinter)/(mu0*Ms);
    Hcrit=0;
    ramping_time=35e-12;
    HSO=10*(10^4/12.54);

    
    if t < ramping_time
        y = max(HSO/ramping_time*t - Hcrit, 0);
    elseif t < 45e-12
        y = HSO-Hcrit;
    else
        y = 0;
    end
    
    dxdt(1) = x(2);
    dxdt(2) = -C1*x(2)+C2*Delta*y;
    dxdt = dxdt(:);
end

As you can see, Delta appears in the definition of dxdt(2). On the other hand, the code which calls this odeFun function is the following:

data=load('data.dat');
ic = [0; 0];
time=data(:,1); % s
Delta=data(:,2).*(10^(-9)); % m
[~,sol] = ode45(@(t,x) odeFun(t,x,Delta),time,ic);

This does not work. It only works if the following lines the definition of the function are rewritten:

function dxdt = odeFun(t, x, Delta, time)
dxdt(2)=-C1*x(2)+C2*interp1(time(:),Delta(:),t)*y;

I do not know if this do exactly what I want. I mean, in principle I do not need to interpolate nothing, because each value of Delta corresponds to each value of the time array. Does this do what I want or there is other way?

Answer 1

[Bjorn Gustavsson]

Your problem is that ode45 will evaluate odeFun at arbitrary points in time (as seen from your perspective at least). That means that odeFun will have to have some way to handle the variable Delta in continuous time - while you only have estimates at some kind of discrete times, that is what you achieve with the interpolation-step.

What you have to determine is how your estimated Delta should be converted to continuous time. This is where it becomes tricky for anyone to give solid detailed advice, since it is not obvious what your ODE represents (it looks vaguely like an equation of motion, but might be something else) nor do we know how you've estimated Delta. If you for example have Delta as some sort of constant force during time-intervalls between your time-steps you could change the interpolation from 'linear' to 'nearest' (you might have to shift the ramping_time by half a time-step), or if you think that Delta should vary more smoothly you might change the interpolation-method to 'pchip' or 'cubic'. (or you could compare all 3-4 interpolation-methods to see how different solutions you get.)

So outside of general advice like this it will be up to you to think about what your Delta-estimate represents - or explain in more details.

HTH

[Roderick]

Hi! Thanks for your answer. Okay, I understand your point. Yeah, it is a equation of motion, where x(1) represents spatial position and x(2) the velocity of an object. Nevertheless, the point is that I didn't want to rely at all in a hypothetical analytical expression for Delta in this approach. What I am thinking about right now is the following. Maybe if I do a spline of my data of Delta to third grade equation, a*t^3+b*t^2+c*t+d, I could extract the parameters a, b, c, d and construct an analytical expression based on that. If that was the case, would it be a legitimate approach or do you think that spurious results could appear in solving my problem because it is such an artificial equation?

[Bjorn Gustavsson]

OK, then Delta is a factor to a force-term. For this it now depends on how it varies with time in your estimate. Will a single cubic polynomial capture all the relevant variation? Is so, by all means use a cubic. If you have more short-time variations then you'll have to use something with more variations. If you want to get a term varying such that it passes through the exact points [time(i_t), Delta(i_t)] then you should use the interpolation. My first lazy stab would be to use the 'pchip' interpolation method:

dxdt(2)=-C1*x(2)+C2*interp1(time(:),Delta(:),t,'pchip')*y;

But if your Delta has some noise-like variability you might prefer some apporximating function, spline, harmonic approximation or something else that capture the physical variation but not the noise - unless the noise is a physical feature.

HTH

[Roderick]

Great! Thank you for your advice! I think that I can continue with those comments.