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

1 次查看（过去 30 天）
Roderick
Roderick 2020-5-22
评论： Roderick 2020-5-22
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?

请先登录，再回答此问题。

采纳的回答

Bjorn Gustavsson
Bjorn Gustavsson 2020-5-22
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
  3 个评论
显示 1更早的评论
Bjorn Gustavsson
Bjorn Gustavsson 2020-5-22
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

请先登录，再进行评论。

更多回答（0 个）

类别

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

标签

Community Treasure Hunt

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

Start Hunting!

Translated by