solution to 2nd order ode having time varying coefficients using RK4 steps in Matlab

Hello everybody, I am trying to solve a 2nd order ODE using RK4. The coefficients are time dependent and assume new values after a fixed time interval and needs to be updated from the solution to previous steps. I think I have to use a loop. Can I connect this to any standard suite like ode45? Or do I need to write down the RK4 algorithm in steps and then use a loop counter?

Note: The fixed time interval can be same as the stepped time interval we use while solving the ode suite.

Generally any IVP evolves with time for a given set of coefficients. In this case if the coefficients themselves were to change every second will I be able to capture the transient dynamics or I need to go for a time step one order lower with the ode suite to capture the dynamics? Then how do I run a loop along with the suite like ode45?

Thank you.

[James Tursa]

Please tell us more about the time varying nature of the coefficients. Are they step functions? Do they happen based on some event? Are they functions of the state vector? Etc.

Can you provide an example of how they change?

[Subrata Chakrabarti]

Dear James,

The problem I am trying to solve is a 2nd order spring mass damper initial value problem with a forcing function. The position, of a valve piston,as determined by the equation in its first step will lead to an area through which flow of gasses take place. This leads to a fall in source pressure which in turn will cause a change in the forcing function for the next step. This becomes the input forcing function for the next step .Hope this explains.

Answer 1

[Ameer Hamza]

ode45 is capable of handling time-varying parameters. See this example for details: https://www.mathworks.com/help/matlab/ref/ode45.html#bu3l43b

You don't need to write a loop yourself. Just modify the parameter inside the ODE function according to your requirement.

[Subrata Chakrabarti]

Thank you @ Ameer Hamza for your response. As a small clarification I would like to add that the time varying coefficients gets updated based on the solution to the immediate previous step. These coefficients are not sine or exponential functions. They do not follow any pattern.

In brief, what I would like to know is how to approach the solution to a 2nd order diff equation using RK4 in matlab with time varying coefficients updated based on previous solution steps.

Thank you once again for your response and look forward to hearing from you soon.

[Ameer Hamza]

Can you clarify what do you mean by "time varying coefficients updated based on previous solution steps"? Can you give an example of the ODE system?

[Subrata Chakrabarti]

Dear Mr Hamza,

Apologies for the dealy in responding.

The ode I am trying to solve is closest in form to the one cited below(quoted from the link you shared).

The only difference is in my case , the function 'f' or 'g' are not interpolated values but they derive their values from previous time steps. They evolves as the solution progresses. Also, the variable 'y' has an upper and lower bound to itself.

ODE with Time-Dependent Terms

View MATLAB Command

Consider the following ODE with time-dependent parameters

The initial condition is . The function f(t) is defined by the n-by-1 vector f evaluated at times ft. The function g(t) is defined by the m-by-1 vector g evaluated at times gt.

Create the vectors f and g.

ft = linspace(0,5,25);

f = ft.^2 - ft - 3;

gt = linspace(1,6,25);

g = 3*sin(gt-0.25);

Write a function named myode that interpolates f and g to obtain the value of the time-dependent terms at the specified time. Save the function in your current folder to run the rest of the example.

The myode function accepts extra input arguments to evaluate the ODE at each time step, but ode45 only uses the first two input arguments t and y.

function dydt = myode(t,y,ft,f,gt,g)

f = interp1(ft,f,t); % Interpolate the data set (ft,f) at time t

g = interp1(gt,g,t); % Interpolate the data set (gt,g) at time t

dydt = -f.*y + g; % Evaluate ODE at time t

Solve the equation over the time interval [1 5] using ode45. Specify the function using a function handle so that ode45 uses only the first two input arguments of myode. Also, loosen the error thresholds using odeset.

tspan = [1 5];

ic = 1;

opts = odeset('RelTol',1e-2,'AbsTol',1e-4);

[t,y] = ode45(@(t,y) myode(t,y,ft,f,gt,g), tspan, ic, opts);

Plot the solution, y, as a function of the time points, t.

plot(t,y)