0

I have a vector x of length N and I want to use its values for solving the differential equation: dy/dt = x - 4*y. For each step I want the ode solver function to use one value of the vector and then use the next value of the matrix for the next step.

I tried doing so by declaring the vector a global variable and use it like this inside the ode solver:

global x
tspan = 0:0.01:10;
[t,filter_coef] = ode45(@ode_filter,tspan,0);

and the solving function like this:

function dx = ode_filter(t,fil)

    global x 
    dx = x - 4*fil(1);

end

Following error was produced

ODE_FILTER returns a vector of length 1002, but the length of initial conditions vector is 1. The vector returned by
ODE_FILTER and the initial conditions vector must have the same number of elements.
Teo Protoulis
  • 203
  • 4
  • 12
  • This question is VERY similar to another one posted very recently: https://stackoverflow.com/q/56206969/3372061. Please see if the solution there works for you. – Dev-iL May 20 '19 at 07:10
  • This is possibly an XY problem (as is often the case when you start hacking interna of standard methods), meaning that you are trying to solve Y for your idea of solving X, so only describe Y. However, your idea of solving X might be sub-optimal, a better solution of X might not require Y or require a modified version of Y. Please tell more of the greater picture of problem X. – Lutz Lehmann May 20 '19 at 08:27

2 Answers2

1

enter image description here

  • y is a function of t, not function of x.

  • From the above equation unless you define x value first, y will be a function of t and x.

  • Also ode45 return only double values, not the function y itself but it's evaluation at different t

  • What can be done is defining a function that throws a given unknown value x to ode45. Then you get the function I deduct above.

  • Next Just set a given x and evaluate the function. The function evaluation gives structure data type which main fields are x and y. x is your t and y remains your y

    global x
    tspan = 0:0.01:10;
    f = @(x)ode45(@ode_filter, tspan,0);
    %% Set x Value First
    x = 2;
    var = f(x);
    t = var.x;
    y = var.y;
Adam
  • 2,726
  • 1
  • 9
  • 22
  • 1
    "_`ode45` excluding `y` and `t` does not accept variable_" this claim is **false**, and I know this by looking at the singnature of `ode45`: `ode45(ode,tspan,y0,options,varargin)`. Notice the `varargin` there? This is MATLAB's way of defining a [variadic function](https://en.wikipedia.org/wiki/Variadic_function), which can accept as many inputs as you want. Please clarify or remove this claim. – Dev-iL May 20 '19 at 07:16
  • Your solution throws an error stating: `ODE_FILTER returns a vector of length 1001, but the length of initial conditions vector is 1. The vector returned by ODE_FILTER and the initial conditions vector must have the same number of elements.` Any way of solving this because as I can see results are calculated perhaps as expected ? – Teo Protoulis May 20 '19 at 16:37
1

You might be basing your approach on your understanding of the explicit Euler method. If that is the case you would be best served by implementing the (exponential) explicit Euler method.

The solver ode45, as all of the other matlab ODE solvers (without specific options?) has variable step size which is automatically adapted to the problem and the current state. Additionally, during each step the passed right side ODE function is evaluated multiple times. Further, the step size regulator depends on the smoothness to a high order of the right side function, non-smooth loci lead to a drastic reduction in step size and possibly multiple restarts from the same current state, further contradicting your assumptions.

Thus, even if you succeed in implementing your idea, you are adding essentially random noise to your function, rendering the solver useless as basic assumptions are violated. Even if a result is produced, it will have almost nothing to do with what you wanted to achieve.

The fastest way to achieve something like what you want to produce is to identify the times at to which the x values are associated and use some interpolation function, zero order hold or linear interpolation, to get the correct value at variable times.

Use solution extension to solve every segment using a smooth right side, changing the constant x(k) for each segment [ tx(k), tx(k+1) ]. Define the function to have a parameter, avoiding the hassle of global variables

function dy = ode_filter(t,y,x)
    dy = x - 4*y;
end

then call the integrator for the first segment to initialize and then for all the remaining segments using their constant and the segment end.

sol = ode45(@(t,y)ode_filter(t,y,x(1)), [ tx(1) tx(2) ], y0)
for k in 2:N
    sol = odextend(sol,@(t,y)ode_filter(t,y,x(k)),tx(k+1));
end

(Always think about using the options mechanism to set costumized, problem specific error tolerances.)

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51
  • Does it change anything if I am producing the x vector from another ode solving function ? – Teo Protoulis May 20 '19 at 08:18
  • The best would be if you would then solve the coupled system as a coupled system. If you use the matlab integrators for the other integration, use the `sol_x=ode45...` object return option, this gives you an automatically method adapted interpolator in `x=deval(sol_x,t)`. – Lutz Lehmann May 20 '19 at 08:22
  • I am using ode45 to solve the other equation: [t,x] = ode45(...). This is how I produce the x vector. – Teo Protoulis May 20 '19 at 08:27
  • Is there a requirement of only using the sampled data? Of using a specific interpolation method? You could use `function dxy = sys(t,xy) dxy = [ f_x(t,xy(1)), f_y(t,xy(1),xy(2)) ]; end` and then `[T,XY] = ode45(@sys,...)` to solve the coupled system in one go. – Lutz Lehmann May 20 '19 at 08:31
  • It would be best to use matlab’s ode integrators. I have implemented the solution using euler’s method for every differential but I have to implement it using the ode integrators. The [f_x(t,xy(1)) , f_y(t,xy(1),xy(2)) ] vector represents the equations I have to solve ? – Teo Protoulis May 20 '19 at 08:35
  • Yes, they stand for the system `dx/dt = f_x(t,x), dy/dt = f_y(t,x,y)` (not partial derivatives, just names), I should have put that as intro in the last comment. – Lutz Lehmann May 20 '19 at 08:56
  • If I use this type of system you suggest, how do I pass the initial condition of each variable ? Using a vector as normally ? – Teo Protoulis May 20 '19 at 09:03
  • Yes, everything is vectorized, `xy0 = [x0, y0]`. – Lutz Lehmann May 20 '19 at 09:18
  • Will try them and come back with feedback. – Teo Protoulis May 20 '19 at 09:19
  • Putting all the differential equations in one system and solving them at one go works but results are not at all what was expected. – Teo Protoulis May 20 '19 at 17:45
  • Then you would need to explain what the reference values are and how they were obtained. As you have the word "filter" in your text, is there some time-frequency transformation involved? Convolution operators related to some discretization scheme? If this leads to massive changes to the question text, open a new question instead. Perhaps the overarching topic needs to move to signal processing (dsp) or scientific computing (scicomp) on stackexchange.com – Lutz Lehmann May 20 '19 at 17:49
  • Task is about implementing steepest descent for online parameter estimation of a first order system. I have to estimate two parameters. The word filter refers to a first order filter which filters the state and input vectors of the system in order to linearly paremeterize it and then apply the steepest descent method to estimate the parameters. Overall I have 5 differential equations: 1 produces the state vector, 2 produce the vectors of the fiend 2 the vectors of the estimations. – Teo Protoulis May 20 '19 at 17:53
  • That looks like something that, if you want to discuss it at all, is better placed on scicomp – Lutz Lehmann May 21 '19 at 06:48