Convergence with fsolve?
19 次查看(过去 30 天)
显示 更早的评论
Hi all,
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.
% Inputs
% Load Functions
load('Struct_Rev00.mat')
% Engine Data
Nc = 9 ; % Number of Cylinders
% Solver Inputs
N = length(J) ; % Degrees of Freedom
Phi = 2 * pi / Nc * (0:2)' ; % DOF Phase angle
% Solve System
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) ;
%sysfcn(thsol, J, K, Tgf, Tm2f, Phi, N)
% System of Equations Function
function Fu = sysfcn(th, J, K, Tgf, Tm2f, Phi, N)
% Preallocate
Fu = zeros(N, 1) ;
% Equation RHS
Theta = wrapToPi(th(1:N-1) - Phi) ; % Angle argument
T_2 = [reshape(Tgf(Theta), N-1, 1) ; 0] ; % T_2 function
T_1 = diag([reshape(Tm2f(Theta), N-1, 1) ; 0]) % T_1 function
% Equation
Fu(1:N) = (J - T_1) \ (-K * th(1:N) + T_2) ;
end
Thanks for your help in advance,
KMT.
0 个评论
回答(0 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!