Hello, I want to solve system of non linear equations using fsolve. There are thre equations and two unknowns K and L. looking at my code, please advice how to best do this. Also please advice if this is the best way or there is an alternate option.

Ld= 0.8194  %constant
% cos^2(K*L)+ Ld*K*sin(K*L)cos(K*L)=2.2 two non linear equations
% cos^2(K*L)=0.299
% Ld*cos^2(K*L)+sin(K*L)cos(K*L)/K=0.262

cos^2(K*L)+ Ld*K*sin(K*L)cos(K*L)-2.2=0
cos^2(K*L)-0.299=0
Ld*cos^2(K*L)+(sin(K*L)cos(K*L))/K-0.262=0

function F=calib(K,L)
F(1)=cos(K*L)^2+ Ld*K*sin(K*L)*cos(K*L)-2.2;
F(2)=cos(K*L)^2-0.299;
F(3)=Ld*cos^2(K*L)+(sin(K*L)cos(K*L))/K-0.262;

fun = @calib (K,L);
initial_val=[2.73,0.6]  % initial value of K and L respectively
x = fsolve(fun,initial_val);

[Star Strider]

In the future, please highlight the code in your questions or comments, then click on the {}Code button to format it.

[Sumera Yamin]

Thank you for the guidance. in future i will take care about coding format

Answer 1

[Star Strider]

You have the correct approach and are close to the correct solution. There are some errors in your code for ‘calib’ that I corrected. (Check to be certain that they do what you want.) I also created it as an anonymous function, for convenience.

The optimization functions all expect a vector of parameters, not separate parameters. I changed ‘fun’ to provide the correct calling convention, without having to change ‘calib’.

This runs:

Ld= 0.8194  %constant
calib = @(K,L) [cos(K*L).^2 + Ld*K*sin(K*L).*cos(K*L)-2.2;  cos(K*L).^2-0.299;  Ld*cos(K*L).^2 + (sin(K*L).*cos(K*L))/K-0.262];

fun = @(b) calib(b(1),b(2));
initial_val=[2.73, 0.6]  % initial value of K and L respectively
[x,fval] = fsolve(fun,initial_val);

You may need to experiment with it to actually solve your equations, since the end values are not near zero.

[Sumera Yamin]

thank you for the correction

[Star Strider]

As always, my pleasure.

[Sumera Yamin]

hi, continuing with the same, i have another question. if we increase the no of non-linear equations to lets say 16 with the same unknowns, this should increase the accuracy of my results. but when i try to do this, the programe gives me an error. please let me know what is the error in my programme.

if true

ld= 0.8194;
calib = @(K,L) [cos(K*L).^2- ld*K*sin(K*L)*cos(K*L) - 2.2;
ld*cos(K*L).^2+sin(K*L)*cos(K*L)/K - 0.262;
sin(K*L)*cos(K*L)-ld*K*sin(K*L).^2 + 1.21;
sin(K*L).^2/K + K*sin(K*L)*cos(K*L) - 0.63;
-K*sin(K*L)*cos(K*L) - 2.15;
cos(K*L).^2 - 0.299; 
-K*sin(K*L).^2 + 2.13;
sin(K*L)*cos(K*L) - 0.7;
ld*K*sin(K*L).^2-sin(K*L)*cos(K*L) - 1.21;
sin(K*L).^2/K-ld*K*sin(K*L)*cos(K*L) + 0.63;
cos(K*L).^2- ld*K*sin(K*L)*cos(K*L)- 2.22;
ld*cos(K*L).^2 + sin(K*L)*cos(K*L)/K - 0.299;
K*sin(K*L).^2- 2.13;
sin(K*L)*cos(K*L)+ 0.7;
-K*sin(K*L)*cos(K*L)- 2.15;
cos(K*L).^2- 0.299]; 

fun = @(b) calib(b(1),b(2));
initial_val=[2.73, 0.06];  % initial value of K and L respectively
[x,fval] = fsolve(fun,initial_val)

if true
  % code
end

[Star Strider]

It does not give an error. It gives a Warning and that it was unable to converge on a solution:

No solution found.

fsolve stopped because the last step was ineffective. However, the vector of function
values is not near zero, as measured by the default value of the function tolerance. 

<stopping criteria details>

x =

     2.654317999304830  -0.323464888160914

and a vector of final values for the function, none of which are close to zero.

There is no error in your code. Your system does not have a solution, at least near your initial parameter estimates.

I would instead use an unconstrained genetic algorithm (the Global Optimization Toolbox ga (link) function) if you want to search the entire parameter space for solutions. Experiment with the various options and initial populations. (I always define my own random initial population of at least 1000 individuals.) There are other global search function options as well.

Note that a global parameter search may not yield a solution, since no solution may exist for your system.

[Sumera Yamin]

ok, but i do not understand what is mathematically wrong with the system that it doesnot have a solution. So in theory, "K" and "L" represents real physical quantities for which we know the estimated initial values. The whole problem goes like this. I have my data x and y, i plotted the data does the curve fitting and the equations above are transfer functions of my data. So, i want to calibrate the "K" and "L" values which fits my data according to above equation. This is explained in the following code. For simplicity, i have again reduced the no of equations to 3. Can you please comment what is mathematically wrong with my system?

if true
  % code
data=[0  0;
  0.05 1.108646244630E-01;
  0.10 2.217423074817E-01;
  0.15 3.325947375398E-01;
  0.20 4.434863433851E-01;
  0.25 5.543595496420E-01;
  0.30 6.652338361973E-01;
  0.35 7.761094191116E-01;
  0.40 8.869865144820E-01;
  0.45 9.978653384221E-01;
  0.50 1.108746107036E+00];

k=10.9251;                   %fixed
xdata=data(:,1);
ydata=data(:,2);

% where data is related by the following euation.
xdata=(cos(K*L).^2- ld*K*sin(K*L)*cos(K*L))*ydata;
%when i plot x and y and do linear curve fit as shown in attached fig, the slope is 2.2. So,
cos(K*L).^2- ld*K*sin(K*L)*cos(K*L) = 2.2;  % 2.2 is the slope from the linear curve fitting to the above data
%Similarly for othersets of data yields following equations
%cos(K*L).^2 = 0.299; 
%Ld*cos(K*L).^2 + (sin(K*L).*cos(K*L))/K = 0.262;

% so i intend to solve these 3 equations using fval to get the Value of K and L from
% data. My sloution should be closed to values given in initial_val. This
% time use initial value as given below.
 initial_val=[5.46,0.14]

end

[Star Strider]

There may not be anything ‘mathematically wrong’ with your system. There may not be a solution, especially since ‘K’ and ‘L’ are arguments to trigonometric functions.

If you create a matrix version of ‘calib’:

calibm = @(K,L) [cos(K.*L).^2 + Ld*K.*sin(K.*L).*cos(K.*L)-2.2*ones(size(K));  cos(K.*L).^2-0.299*ones(size(K));  Ld*cos(K.*L).^2 + (sin(K.*L).*cos(K.*L))./K-0.262*ones(size(K))];

and then plot contours of each component matrix:

Nv = 50;
v = linspace(0,1,Nv);
[Km,Lm] = ndgrid(v*10, v);

Z = calibm(Km,Lm);

figure
contour(Km, Lm, Z(1:Nv,:), [0 0], '-r')
hold on
contour(Km, Lm, Z((1:Nv)+Nv,:), [0 0], '-g')
contour(Km, Lm, Z((1:Nv)+2*Nv,:), [0 0], '-b')
hold off
grid on
xlabel('K')
ylabel('L')
legend('F(1)', 'F(2)', 'F(3)', 'Location','SW')

you will see that at no one point do all of the zero-contours have the same ‘K’ and ‘L’ values. Your system simply has no solution, at least in this region. If you do a global search, you may be able to find a solution. I leave that to you to experiment with.

(The code for the plots is not efficient. However it is a ‘one-off’ designed to demonstrate the problem, and is not necessary for the optimization.)

[Sumera Yamin]

alright, code for the plot is not working for me as my MATLAB is giving me some license issue for the symbolic toolbox. I will now check for the global search. Thank you very much for the explanation.

[Star Strider]

My pleasure.

My code does not use the Symbolic Math Toolbox. You should be able to run it without problems.

[Sumera Yamin]

yeah, so i did got the contours but i see, the red line did cross other blue and green for two values? how can i explain? I still cannot digest this concept because, theoretically my physical system should have some real solution.

[Sumera Yamin]

I just noticed the red line crossing the green and blue for K between 5-6 and L=0.1 and 0.2. this is the range of K and L which i am expecting as a solution.

[Star Strider]

The lines cross each other, however they do not all cross at the same points.

If you want to use the ga or other optimization functions, your objective function will have to return a scalar value. One way to so this is to use the norm function with ‘fun’:

fun = @(b) norm(calib(b(1),b(2)));

I cannot claim that is mathematically correct for the analysis you want to do. I can claim that it will work with ga and other functions that require that the objective function to return a scalar result.

[Sumera Yamin]

sorry for the stupid question, but what is ga?

[Star Strider]

Not stupid at all. The ga (link) function is the genetic algorithm implementation in the Global Optimization Toolbox.

[Walter Roberson]

You appear to be trying to match experimental data to theoretical equations as if each entry were exact and the model was perfect. You should instead be looking for parameters that minimize the overall error. For that task, I recommend the Curve Fitting Toolbox.