ODE solver, how to pass a vector value used in the differential equations

Hi,

I am solving a system of differential equations of the following format:

dy = single(t, y, J)

dy(1) = (A*(y(3))/(B*y(1)))*y(1)+ C*y(3)^2; 
dy(2) = (B*(y(3))/(A*y(1))-1); 
dy(3) = J/C*y(3) - D*(y(3))/(A*y(1))*y(1);

where A,B,C,D are all constants.

and I am initiating the solver as

[t_ ,y_]=ode45(@single,timespan,y0,options, J);

I am defining my initial condition all as zeros in vector y0, and my timespan is a vector defining the times at which I require my results. If I designate a constant for J, I am getting the result I am expecting.

but in reality, "J" should be a time_varying variable which has same size as "timespan", and is placed in a vector, thus I require my ODE45 to solve my system at timespan(x) where the value of J is J(x). Thus I just wanted to know how I can acheive this.

Regards, Arsalan

[MOHAMMADAMIR KARIMZADEH]

Hi Arsalan, I am having a small problem while solving a system of ODE similar to the one you mentioned above.cod you please share your matlab code for to solve the above system of ODE.

[Walter Roberson]

See my just-posted response https://www.mathworks.com/matlabcentral/answers/52740-ode-solver-how-to-pass-a-vector-value-used-in-the-differential-equations#answer_704728

Answer 1

[Jan]

Shadowing the built-in function single() can have strange effects. Better use a different name.

nTime     = length(timespan);
y_list    = zeros(nTime, length(y0));
y_list(1, :) = y0;
for iTime = 2:nTime
  [t_, y_] = ode45(@yourFcn, timespan(iTime-1:iTime), y0, options, J(iTime));
  y0       = y_(end, :);
  y_list(iTime, :) = y0;
end

Now y_list contains the y-Value for the time in timespan.

Another approach would be to check inside the function to be integrated the value of t and use a corresponding J value. But this would result in discontinuities and the stepsize control of ODE45 would either explode or lead at least to results, which are dominated by rounding errors. Although there are some systems, which can be integrated "successfully" in spite of discontinuities, it is a bad idea to drive numerical software apart from the specifications. Therefore the piecewise integration as shown above is the only scientifically correct solution - some mathematicians even claim, that numerics cannot be counted to science at all...

[Ole]

This is a nice implementation but is extremely slow can it be vectorized ?

[Jan]

@Ole: No, this cannot be vectorized or parallelized, because the initial value of a time step depends on the final value of the previous step.

You can use the profiler to identify the bottlenecks. Perhaps your system is stiff and ODE45 is the wrong integrator.

[DIP]

Sir, how do i solve an ode with an array , using ode45. ODE is udY/dx=P, Y,P are arrays. is it the same approach above ?

[Teo Protoulis]

The solution I get from solving the differential equation is placed inside y_list and not y_ ?

[mohammad heydari]

hi jan.

Assuming j have a length of 2001 and tspan has a length of 50001, what would you suggest?

I have such a problem in my schedule.

Regards, mohammad

Answer 2

[Arsalan]

Thanks Jan, It works like a charm. your right about the the claim that "numeric cannot be counted to science" but in this case we have no other choice since dt is not equal to 0. so we just have to accept

Answer 3

[Walter Roberson]

If you have a variable (such as J in this example) that is defined at particular time steps, then mathematically each them represents a discontinuity to the equations (or at least to the derivatives of the equations), and that is a problem for ode*() solvers.

If you try to do linear interpolation of the variable according to the current time, then that will still have discontinuity of derivatives.

If, however, you can do spline interpolation of the variable, then that has the needed mathematical properties.

One way to do spline interpolation of a constant vector efficiently is to use mkpp() to construct a piecewise polynomial, and pass the coefficients of that into the function and inside the function, use ppval() to do the interpolation.

https://www.mathworks.com/help/matlab/ref/mkpp.html

A couple of notes about this approach, though:

1. This approach is completely unsuitable for cases where you have an instantaneous (or near instantaneous) impulse. It is only suitable for cases where the variable effectively represents a function that has been sampled at locations and is to be interpolated at other locations
2. This kind of interpolation is done forwards and backwards in time. If the data is [1, 4, 3] at t = [0, 1, 2], then at t = 0.99 (almost 1) then the interpolated value would be something between 1 and 4 that was almost 4 -- so the value at t = 1 is "smeared" back in time and forward in time both. If you need a stepwise sequence, such as the data is to be strictly 1 between t = 0 and t = (1 minus epsilon) and then is to suddenly become 4 at t = 1 exactly, then this kind of approach cannot be used.
3. For example this approach is not at all usable for cases where you are modeling hitting a surface and bouncing.

In all other cases (instantaneous impulses for example) then there are two ways to proceed:

• if the changes are at specific times, then call the ode*() function giving a tspan that ends at the time of the new value. That will stop the integration there. Then make whatever changes are appropriate to the boundary conditions (such as adding the instant impulse), and then call the ode*() function with the modified boundary conditions. Modifying the boundary conditions might not be required. If for example a heater switches on after 10 seconds, then tspan [0 10] the first time, and take the resulting conditions as initial conditions for anoher ode*() . The general principle is that it is okay to use conditional statements inside your ode function provided that only one branch is used in any one ode*() call (or if you have been very careful to match the derivatives at the boundary conditions.)
• If the changes are not at specific times, then call the ode*() function passing in an options structure that defines an event function that notices the condition and terminates the call; then outside of ode*(), modify the boundary conditions (if required) and restart from the current time. See the ballode example.