Index exceeds matrix dimensions, using fsolve

I don't see any reason why I should be exceeding matrix dimensions here. The constants in the function (z onwards till Dinf) are defined in a separate script,which I 'run' just before executing the function.

function [solution]=nano2(x,z,PHI,gamma,Kc,Dp,PHI_B,Mc,N,N_c,deltax_j,J_v,f,gc,T,Dinf)
A=x(1:3,1:N+2); %Number of rows to be modified with number of components
B=x(4:6,1:N+2);
Q=x(7:9,1:N+2);
R=x(10:12,1:N+2);
E=x(13:15,1:N+2);
F=x(16:18,1:N+3); %Might need to look at its dimensions.
PSI=x(19,1:N+3);
c=x(20:22,1:N+3);
zeta=x(23,1);
for c1=1:N_c
for c2=3:N+2
solution=[A(c1,c2)-((Dp(c1)/deltax_j)-(1/2)*Kc(c1)*J_v+(1/2)*z(c1)*Dp(c1)*(f/(gc*T))*((PSI(c2+1)-PSI(c2))/deltax_j));%Membrane active layer domain for ion i-linearization of Nernst-Planck equation--------(25-31)
B(c1,c2)+(Dp(c1)/deltax_j)-(1/2)*Kc(c1)*J_v+(1/2)*z(c1)*Dp(c1)*(f/(gc*T))*((PSI(c2+1)-PSI(c2))/deltax_j);
Q(c1,c2)-(1/2)*z(c1)*Dp(c1)*(f/(gc*T))*((c(c1,c2+1)+c(c1,c2))/deltax_j);
R(c1,c2)+(1/2)*z(c1)*Dp(c1)*((f/(gc*T)))*((c(c1,c2+1)+c(c1,c2))/deltax_j);
E(c1,c2)+J_v;
F(c1,c2)+(1/2)*z(c1)*Dp(c1)*(f/(gc*T))*(c(c1,c2)+c(c1,c2+1))*((PSI(c2+1)-PSI(c2))/deltax_j);
F(c1,c2)-A(c1,c2)*c(c1,c2)-B(c1,c2)*c(c1,c2+1)-Q(c1,c2)*PSI(c2)-R(c1,c2)*PSI(c2+1)-E(c1,c2)*c(c1,N+3);
% Feed/solution equilibrium
%Feed solution/membrane mass transfer coefficient for each component------(22)%
(Mc(c1)-J_v)*c(c1,2)+J_v*c(c1,N+3)+(z(c1)*c(c1,2)*Dinf(c1)*(f/(gc*T)))*zeta-Mc(c1)*c(c1,1);
%Thermodynamic equilibrium conditions at feed solution/membrane for each component--------(24)
(gamma(c1,3)*c(c1,3))+(gamma(c1,2)*c(c1,2)*z(c1)*(f/(gc*T))*PHI(c1)*PHI_B(c1)*exp((-z(c1)*f*PSI(3))/(gc*T))*PSI(3))-(gamma(c1,2)*c(c1,2)*PHI(c1)*PHI_B(c1)*exp((-z(c1)*f*PSI(3))/(gc*T))*(1+(((z(c1)*f*PSI(3))/(gc*T)))));
%Thermodynamic equilibrium conditions at membrane/permeate-solution interface for each component----(33)
(gamma(c1,N+2)*c(c1,N+2))-(PHI(c1)*exp(z(c1)*(f/(gc*T))*(PSI(N+3)-PSI(N+2))))*c(c1,N+3)+(gamma(c1,N+3)*c(c1,N+3)*z(c1)*(f/(gc*T))*PHI(c1)*exp(z(c1)*(f/(gc*T))*(PSI(N+3)-PSI(N+2))))*PSI(N+2)-(gamma(c1,N+3)*c(c1,N+3)*z(c1)*(f/(gc*T))*PHI(c1)*exp(z(c1)*(f/(gc*T))*(PSI(N+3)-PSI(N+2))))*PSI(N+3)-(gamma(c1,N+3)*c(c1,N+3)*PHI(c1)*exp((z(c1)*(f/(gc*T)))*(PSI(N+3)-PSI(N+2))))*(z(c1)*(f/(gc*T))*(PSI(N+3)-PSI(N+2)));
%Electroneutrality at each node--------(32)"
z(1)*c(1,c2)+z(2)*c(2,c2)+c_x;
%Electroneutrality in boundary layer on feed solution/membrane------(23)
z(1)*c(1,2)+z(2)*c(2,2)+z(3)*c(3,2);
%Electroneutrality in boundary layer on permeate solution/membrane------(34)
z(1)*c(1,N+3)+z(2)*c(2,N+3)+z(3)*c(3,N+3);
];
end
end

end

Answer 1

[Star Strider]

I suggest you not worry about fsolve just now. Call your nano2 function directly with your initial parameter estimates for it, then see if you can isolate the statement that is likely throwing the error. Also, look at the sizes of A through zeta and check c1 and c2.

The debugger will likely be a significant help in code as complex as this.

[Yagnaseni Roy]

Thanks for the advice, Star Strider. However, in waiting for an answer to this question, I tried out taking out one line each from inside the fsolve equation set (as shown above) and I was able to detect the error. MATLAB can also be unpredictable, because I made the same changes in another similar file and it gave errors. Anyway, I think that for now I have this one sorted.

However, could you tell me if it is possible to do iterations on fsolve? Even though I got answers, they don't look correct. I feel like it might be coz I need more iterations. Is it possible to do so in fsolve?

Thanks for the help!

[John D'Errico]

More iterations will simply not do anything.

If fsolve has converged, letting it do more iterations is a waste of effort, as it will stop in the same place.

System of equations often have multiple solutions, and you may have converged to the wrong one. Better starting values are crucial here. After all, if you don't like the answer, then you probably know what the parameters should look like.

Oh, and DEFINITELY learn to use the debugger. It would have saved you a great deal of time and effort.

[Star Strider]

My pleasure!

You can change fsolve’s default behaviour with optimoptions. To increase its iterations with fsolve, I’d do something like:

opts = optimoptions('fsolve', 'MaxFunEvals',10000,  'MaxIter',10000, 'TolFun',1E-8);

Then ‘opts’ becomes the third argument in your fsolve call. It should converge before it reaches those limits. Also, experiment with several different starting parameter estimates, even extreme ones on the order of 1E±5. (Unless you have reason to believe they are larger than that, I would suggest those as the maximum magnitudes to avoid losing precision with powers of parameters.) I have no idea what your regression hypersurface ‘looks’ like, but there are likely many local minima. You may have to wander a bit to find the correct one.

EDIT — John is correct, however I assumed from your comment that fsolve had stopped at its default iteration limits before it converged.

[John D'Errico]

It may matter if fsolve has stopped having exceeded the maximum number of iterations. If it is just your vague feeling, then don't bother. Giving it more rope then will not help a bit.

So, IF fsolve has stopped at the maximum iterations, then it MAY help. But the fact is, very often it will not do so at all. When solvers get into a problem like this, they are often in the process of getting hung up in a morass, with no escape. (Maybe not a very good metaphor, but I think it gets the point across.) What I'm trying to say here is most of the time when you get that message and you give it more rope, it will still probably just fuss around and leave you at the same place. It really is a rare thing that more iterations will help fsolve to magically escape from the bad place.

So feel free to try upping the max iteration limit, but just don't waste a lot of time there. I would very much recommend trying to find better starting values.

If all such sets of starting values end up in the same bad place, THEN you need to consider if there is a problem, a reason why this is going to a bad place. I have no idea what this system represents, but consider if there may be a coding problem. Consider if there may be a problem in whatever system this is supposed to model. And finally, consider if perhaps the solution may just be as good as you can expect.

[Star Strider]

I would agree, but Yagnaseni Roy is fitting about 45 parameters in this model. That may require giving the solver more latitude, and time.

[Yagnaseni Roy]

Thanks both of you! I will look into the various options you mentioned and see what works for me.