Changing constant with timespan in ode45 solver

I am writing a code that solves a set of equations from ode45. However, I have a constant that changes in the equations at every time point. E.g.

dA = J(1)*x;

dB = J(5)*x;

So at t = 0, x =0 but at t = 1, x =0.25 etc. The solutions are then used as initial values for the next ode solver loop etc. How can I implement this so the ode can solve with this in mind?

I have tried for loops but still no luck.

Answer 1

[Alan Stevens]

Make x a function of t and call it from the function defining the rate equations ode45 is calling..

[Koren Murphy]

Hi Alan

The variable X does not actually increase at a linear rate - it is defined from a set of interpolated results. It is a vector of numbers (see below - these are simplified values for ease) that do not have a defined relationship. Can you explain further what you mean if it would still apply to this? Thank you!

0

-0.0007

0.0017

0.0043

0.0069

0.0095

0.0121

[Alan Stevens]

Yes, if I understand you correctly, you can use an interpolation function to find X as a function of t (look at interp1 or spline).

[Koren Murphy]

No I mean that is how i have found X. X is a vector of random numbers and I need to use different one in the entire set of equations at each time point.

[Alan Stevens]

In that case put your call ode45 (or whatever one you are using) within a for loop and pass a different value of X to the rate equations each time through the loop.

[Alan Stevens]

Structured something like the following perhaps

% X = [x1 x2 x3 ...etc];
% tspan = [0 tfinal];
% IC = {A0 B0];
% for i = 1:numel(X)
% [t, AB] = ode45(@rate,tspan,IC,[],X(i));
% ...
%  function dABdt = rate(t,AB,x)
%            A = AB(1);  % I'm assuming the J's involve A and B somehow
%            B = AB(2);
%          .....
%          dABdt = [ J(1)*x;
%                    J(5)*x];
% end

[Goutham Sajja]

Dear Alan Stevens,

I also have a similar problem, wherein I am in the process of evaluating the value of inductance of a solenoid for each plunger position, according to the formula:

f = 0.5 * i^2(dL/dx)

i is the current supplied to the solenoid, f is the force on the plunger at each step of its position.

I need to extract the value of L as a function of x.

So, I tried a Matlab code something similar to what you have suggested. Could you please give me a your inputs?

Thanks and Regards,

Goutham Sajja

L0 = 3.8184e-3;                         % inital value of induactance
i = 8.75;                               % current in ampere
x = plungerPosition;                    % data exported in .mat file from excel
for j = 1:numel(F)                      % F is the force at each plunger position, also exported in .mat file from excel
[x,L] = ode45(@eval, x, L0,[],F(j));
end

function dLdx = eval(x,F,i)
dLdx = (F/0.5*i^2);                     % equation is dL/dx = F/0.5*i^2
end
% % 

[Alan Stevens]

If plungerPosition contains a vector of x-values, and you want to pass i to the function, then you shoud probaby replace

for j = 1:numel(F)                      % F is the force at each plunger position, also exported in .mat file from excel
[x,L] = ode45(@eval, x, L0,[],F(j));
end

function dLdx = eval(x,F,i)
dLdx = (F/0.5*i^2);                     % equation is dL/dx = F/0.5*i^2
end

with something like

L0 = 3.8184e-3;                         % inital value of induactance
i = 8.75;                               % current in ampere
x = plungerPosition;                    % data exported in .mat file from excel
for j = 1:numel(F)                      % F is the force at each plunger position, also exported in .mat file from excel
[x,L(:,j)] = ode45(@(x,L) eval(x,L,F(j),i), x, L0);
end
function dLdx = eval(x,L,F,i)
dLdx = (F/0.5*i^2);                     % equation is dL/dx = F/0.5*i^2
end

This is off the top of my head! I can't test it, because I don't have the x and F data.

[Goutham Sajja]

Dear Alan,

I thank you for your reply. With the modified code that you suggested, I get a square matrix for L rather that a column matrix.

I could provide you with the x and F data.

x (mili meter) - 3.5 3.0 2.5 2.0 1.5 1.0 0.5 0

F (newton) - 34.7 44.7 59.2 82.5 124 203.3 331 507.2

Could you please try and give me your inputs?

Thanks and Regards,

Goutham Sajja

[Alan Stevens]

Ah, you mean like this:

L0 = 3.8184e-3;                         % inital value of induactance
xspan = [0 3.5];
[X,L] = ode45(@eval, xspan, L0);

plot(X,L),grid
xlabel('x'),ylabel('L')

function dLdx = eval(x,L)
i = 8.75;                               % current in ampere
X = [3.5 3.0 2.5 2.0 1.5 1.0 0.5 0];
F = [34.7 44.7 59.2 82.5 124 203.3 331 507.2];
f = interp1(X,F,x);
dLdx = f/0.5*i^2;                     % equation is dL/dx = F/0.5*i^2
end

You could use a different interpolation function if interp1 isn't accurate enough for you.

[Torsten]

dLdx = f/(0.5*i^2)

[Goutham Sajja]

Dear Alan,

Thank you so much. This seems to work and give me the value of L as a function of x.

Let me tell you I only gave a sample data of 8 values of plunger position and the corresponding force values.

I have the physical measured data of x and F. The size of each column matrix is 984x1.

In that case, it is not necessary to use the interp1 function right?

Thanks

Goutham Sajja

[Alan Stevens]

You will probably still need the interpolation function as the ode solver might well use values of x that aren't in your list in its internal calculations.

[Bhanu Pratap Akherya]

alan stevens,

i have a set of values for the force term "F" (in the equation my"+cy'+ky=F) saved in an excel file which i have recalled using: F=xlsread('l&d.xlsx','O1:O300');

The problem is that i want to recall 1 value of F at different time steps simultaneously but i dont know how to do that, can you guide me? i have used the ODE function as:

[tsol ysol]=ode15s('beam_function',[1:dt:T],y0);

y0=the initial condition, 'beam_function'=the function file i wrote

[Alan Stevens]

You should make this a completely separate thread, in which you also upload the coding you have implemented so far.

[Bhanu Pratap Akherya]

Thanks for the reply alan, I have posted the question can you please help? Here is the link to the question: https://in.mathworks.com/matlabcentral/answers/1436182-changing-values-of-rhs-with-each-time-step-in-ode

Answer 2

[Ameer Hamza]

See this example on the documentation page of ode45(): https://www.mathworks.com/help/matlab/ref/ode45.html#bu3l43b. It shows how to deal with time-varying parameters.