Error in ode45 (Index exceeds matrix dimensions)

Can someone please explain what I have to fix for the AF_function_in(t) of my code?

My function is as followed:

    function dx=myeqn(t,x,p)
    global tu;

    function output=AF_function_in(t)
      output=interp1(tu(:,1),tu(:,2),t);
      end
    dx(1,1)= (p(1)/p(2))*(AF_function_in(t)-x(1))- (p(4)/p(2))*(x(1)-x(2));
    dx(2,1)=  (p(4)/p(3)) *(x(1)-x(2));
    end

When I run:

tv=1x225 matrix;
[t,x]=ode45(@(t,x) myeqn (t,x,[0.1 0.2 0.3 0.4],tv,[0 0]));

I get an error of:

Error in myeqn/AF_function_in (line 4)
        output=interp1(tu(:,1),tu(:,2),t);

Error in myeqn (line 6)
dx(1,1)= (p(1)/p(2))*(AF_function_in(t)-x(1))- (p(4)/p(2))*(x(1)-x(2));

Error in @(t,x)myeqn(t,x,[0.1,0.2,0.3,0.4])

Answer 1

[Star Strider]

First, ‘tu’ has to be an (Nx2) matrix.

Second, do not use global variables! They make your code much more difficult to troubleshoot. Pass ‘p’ and ‘tu’ as extra parameters instead.

Try this:

function dx=myeqn(t,x,p,tu)
output=interp1(tu(:,1),tu(:,2),t);
dx(1,1)= (p(1)/p(2))*(AF_function_in(t)-x(1))- (p(4)/p(2))*(x(1)-x(2));
dx(2,1)=  (p(4)/p(3)) *(x(1)-x(2));
end

and then call it in your ODE solver as:

[t,x] = ode45(@(t,x) myeqn(t,x,p,tu), tspan, x0);

with ‘x0’ being your (2x1) initial conditions vector.

[Arbol]

Thank you, I have been working too long on this and i get stuck ahha.

and my final code is as followed:

function dx=myeqn(t,x,p,tu)
  output=interp1(tu(:,1),tu(:,2),t);
  dx(1,1)= (p(1)/p(2))*(output-x(1))- (p(4)/p(2))*(x(1)-x(2));
  dx(2,1)=  (p(4)/p(3)) *(x(1)-x(2));
  end

[Star Strider]

My pleasure.

We have all had that experience!

[Arbol]

@Star Strider Also, I see that you understand how to use lsqnonlin or lsqcurvefit by following this: https://www.mathworks.com/matlabcentral/answers/254566-how-do-i-fit-coupled-differential-equations-to-experimental-data.

However, I couldn't able to get my code work. Hopefully you can enlighten me. I have tried all possible way. Using this code I made, the function stopped prematurely, and the First-order optimality converge to one value. I use lsqnonlin function to fit my model as followed:

Fp= 0.1; fp=0.1; fis=0.1; PS=0.1;
param0 = [Fp,fp,fis,PS];
fun=@(p) odefnc_1(p,tu(:,1),tu(:,3),tu);
%Then call the lsqnonlin using:
[newparam,resnorm,RESIDUAL,EXITFLAG,OUTPUT,LAMBDA,Jacobian]=...
    lsqnonlin(fun,param0,lb,ub,options);

function diff= odefnc_1(p,texp,x,tu)
F=p(1);
fp=p(2);
fis=p(3);
PS=p(4);
init = [0 0];

[tv,y]=ode45(@myeqn, texp, init);
P=y(:,1);
I=y(:,2);

diff = (I+P)-x;
function dx=myeqn(t,x)
output=interp1(tu(:,1),tu(:,2),t); 
dx(1,1)= (p(1)/p(2))*(output-x(1))- (p(4)/p(2))*(x(1)-x(2));
dx(2,1)=  (p(4)/p(3)) *(x(1)-x(2));
end
end

[Star Strider]

If lsqcurvefit converged to a local minimum (non-global solution), your code works. It simply did not converge on the optimal (global) minimum. If you know about what the optimal parameter set should be, use those as your starting values.

The Global Optimization Toolbox has functions such as patternsearch that can search a wider area of your parameter space, and may be able to help you find a global optimum. Another option is a genetic algorithm (the ga function). These are not fast, but they will often provide a global optimum if one exists.

Otherwise, keep experimenting with different initial parameter estimates until you get one that are close enough so that lsqcurvefit finds your global minimum.

[Arbol]

it basically returns the initial parameter that put in. For example, if I have x0=[1 1 1 1], it will return [1 1 1 1]. I don't think it is going anywhere.

[Star Strider]

I cannot determine what the problem is. Your code appears to be correct.

One possibility is to vectorize these lines, as:

dx(1,1) = (p(1)./p(2)) .* (output-x(1)) - (p(4)./p(2)).*(x(1)-x(2));
dx(2,1) = (p(4)./p(3)) .* (x(1)-x(2));

See the documentation on Array vs. Matrix Operations (link) for a full description.

[Arbol]

I have tried that, but no luck... I even tried lsqnonlin function, or fit function. Or played around with the options of Tolfun and TolX But no luck on my end. When I run examples of other people's work, their seem to work though @.@.

[Walter Roberson]

I had to invent some numbers to run. I do not see the behaviour you describe:

tu = sort( rand(225,3) );
lb = zeros(1,4); ub = 10000 * ones(1,4);
options = [];

Fp= 0.1; fp=0.1; fis=0.1; PS=0.1;
param0 = [Fp,fp,fis,PS];
fun=@(p) odefnc_1(p,tu(:,1),tu(:,3),tu);
%Then call the lsqnonlin using:
[newparam,resnorm,RESIDUAL,EXITFLAG,OUTPUT,LAMBDA,Jacobian]=...
    lsqnonlin(fun,param0,lb,ub,options);

which I used together with

function diff= odefnc_1(p,texp,x,tu)
F=p(1);
fp=p(2);
fis=p(3);
PS=p(4);
init = [0 0];

[tv,y]=ode45(@myeqn, texp, init);
P=y(:,1);
I=y(:,2);

diff = (I+P)-x;
function dx=myeqn(t,x)
output=interp1(tu(:,1),tu(:,2),t); 
dx(1,1)= (p(1)/p(2))*(output-x(1))- (p(4)/p(2))*(x(1)-x(2));
dx(2,1)=  (p(4)/p(3)) *(x(1)-x(2));
end
end

I notice, by the way, that pass four variables into odefnc_1, but you never use the last of them, tu

[Arbol]

I have the same behavior when I run your code. I reedited my code as your request. As followed:

tu = sort( rand(225,3) );
lb = zeros(1,4); ub = 10000 * ones(1,4);
options = [];
Fp= 0.1; fp=0.1; fis=0.1; PS=0.1;
param0 = [Fp,fp,fis,PS];
fun=@(p) odefnc_test(p,tu(:,1),tu(:,3),tu);
%Then call the lsqnonlin using:
 [newparam,resnorm,RESIDUAL,EXITFLAG,OUTPUT,LAMBDA,Jacobian]=...
     lsqnonlin(fun,param0,lb,ub,options);

[tv,y]=ode45(@(t,x) myeqn(t,x,newparam,tu), tu(:,1), [0 0]);
P=y(:,1);
I=y(:,2);
plot(tu(:,1),tu(:,3),'o', tv,I+P);

I used it with:

function diff= odefnc_test(p,texp,x,tu)
F=p(1);
fp=p(2);
fis=p(3);
PS=p(4);
init = [0 0];
[tv,y]=ode45(@(t,x) myeqn(t,x,p,tu), texp, init);
P=y(:,1);
I=y(:,2);
diff = (I+P)-x;
function dx=myeqn(t,x,p,tu)
output=interp1(tu(:,1),tu(:,2),t); 
dx(1,1)= (p(1)/p(2))*(output-x(1))- (p(4)/p(2))*(x(1)-x(2));
dx(2,1)=  (p(4)/p(3)) *(x(1)-x(2));
end
end

[Arbol]

I get this picture, it varies every time I run it. And getting the same 'stopping criteria details'. Basically this is due to the random numbers that are generated. But the fit is not always accurate.

[Arbol]

Do you think it is because of my ODE? or my parameter?

[Star Strider]

I just now noticed that you are calling the lsqnonlin function, not lsqcurvefit. The lsqnonlin function is appropriate here. The problem is that I do not see where you are comparing the results of your objective function to your dependent variable. (If you are, and I am not seeing it because I do not know what your dependent variable is, please post that function, and your call to it, specifically.)

If ‘y’ is your dependent variable, I would anticipate that your cost function (that you pass to lsqnonlin to optimise) would be something like this (for lsqnonlin):

fun=@(p) y - odefnc_test(p,tu(:,1),tu(:,3),tu);

Try something like that to see if you get the result you want.

Otherwise, it might be easier to use lsqcurvefit instead of lsqnonlin, since it incorporates the cost function in its code. As I mentioned previously, if your regression hypersurface has many local minima, lsqcurvefit may have a very difficult task in getting the ‘best’ fit global minimum. If that emerges as a problem, I would consider using a genetic algorithm (the ga function, although simple but effective ones are not difficult to write yourself) or the patternsearch function.

[Arbol]

Sorry, I kind of messed up with the coding there. I rewrote my function for lsqnonlin as:

function diff= odefnc_1(p,texp,tu)
init = [0 0];
[tv,y]=ode45(@myeqn, texp, init);
P=y(:,1);
I=y(:,2);
diff = (I+P);
function dx=myeqn(t,x)
output=interp1(tu(:,1),tu(:,2),t); 
dx(1,1)= (p(1)/p(2))*(output-x(1))- (p(4)/p(2))*(x(1)-x(2));
dx(2,1)=  (p(4)/p(3)) *(x(1)-x(2));
end
end

And ran:

fun=@(p) odefnc_1(p,tu(:,1),tu)-tissData;
% %Then call the lsqnonlin using:
[newparam,resnorm,RESIDUAL,EXITFLAG,OUTPUT,LAMBDA,Jacobian]=...
    lsqnonlin(fun,param0,lb,ub,options);

And got the same result, but 'newparam' is different than 'param0'. I also tried this for lsqcurvefit, and the same approach. LOL. Bad luck for me. Just one comment, my odefunction will return this attached graph

Where tissue is the result of 'param0', and Data is from the data. Do you think the big gap between the Data and Tissue is what cause this?

[Star Strider]

I am not certain what I am seeing in the plot. If the parameters are now changing, that means the code is working as you want it to. The remaining issue seems to be a scaling problem. and that is a property of your differential equations. The gradient-descent algorithms may be getting ‘stuck’ in a very flat area of your regression space. Different initial parameter estimates could help. You will likely get a better fit with ga or patternsearch in this situation.

[Arbol]

Ok, I will try the other fits. Thank you for the guidance.

[Star Strider]

My pleasure.

[Arbol]

Hey @Star Strider, I got it to fit! I just had to pick out a good starting point (param0) for the fit to work! So that is something to note for. Have to pick out a good initial parameter!

[Star Strider]

Congratulations!

I was confident you could do it!