I was trying to solve a system of equations using fsolve that I havent' been succesfull in doing so. So I have the following 4x4 matrix equation:
where is a diagonal positive matrix (so taking its invrese is no problem), and with are numerical sinusoidal functions with period of 120deg.
In attempting to solve for the values of θ I use the function fsolve, and I have set the number of function evaluations at a sufficiently large number so this doesn't cause the algorithm to stop prematurely. After that, at the end of the solver's attempt I get the following message which obviously means that there hasn't been any convergence to a root since the residuals are quite large, yet the steps taking are extremely small:
fsolve stopped because the relative norm of the current step, 6.019541e-13, is less than
max(options.StepTolerance^2,eps) = 1.000000e-12. However, the sum of squared function
values, r = 1.748432e+00, exceeds sqrt(options.FunctionTolerance) = 1.000000e-03.
What I have tried is:
- Changing through all the solvers, but the answer did not improve (levenberg-marquardt seemed to be the fastest however).
- Setting as a constant, but similarly there was no improvement
- Running the solver through a loop with the final answer of the previous step acting as the intial condition for the next.
- Changed the tolerances so that the algorithm doesn't take too small steps (however it still can't solve as r is always greater than sqrt(options.FunctionTolerance).
Below is my code with the .mat file attached that is needed to run the code.
Nc = 9 ;
N = length(J) ;
Phi = 2 * pi / Nc * (0:2)' ;
options = optimoptions('fsolve','Display','iter-detailed',...
'Algorithm', 'levenberg-marquardt', ...
'MaxIterations', 1000, 'MaxFunEvals', 5000,...
'FunctionTolerance', 1e-3,'OptimalityTolerance', 1e-3,'StepTolerance', 1e-3) ;
th0 = zeros(N, 1) ;
[thsol, ~, ~,Jacobian] = fsolve(@(th) sysfcn(th, J, K, Tgf, Tm2f, Phi, N), th0, options) ;
function Fu = sysfcn(th, J, K, Tgf, Tm2f, Phi, N)
Fu = zeros(N, 1) ;
Theta = wrapToPi(th(1:N-1) - Phi) ;
T_2 = [reshape(Tgf(Theta), N-1, 1) ; 0] ;
T_1 = diag([reshape(Tm2f(Theta), N-1, 1) ; 0])
Fu(1:N) = (J - T_1) \ (-K * th(1:N) + T_2) ;
Thanks for your help in advance,