Trouble using fsolve with two non-linear equations and 3 parameter

I am trying to solve for a and d in the equations

• K_H2O = (((2*a)/(6-a))/(((4-2*a-d)/(6-a))*((4-2*a-d)/(6-a))^0.5)*(p^0.5))
• K_OH = (((2*d)/(6-a))/(((2-a-d)/(6-a))^(0.5)*((4-2*a-d)/(6-a))^(0.5)))

Using fsolve. I also need to pass parameters K_H2O, K_OH and p to the function. So far I have

K_H2O = 20.999;
K_OH = 1.0478;
p = 50;
guess = [2;2];
[a,d] = fsolve(eql_coef(K_H2O, K_OH, p), guess)

and in my function "eql_coef.m" I have

function F = eql_coef(K_H2O, K_OH, p)
F(1) = (((2*a)/(6-a))/(((4-2*a-d)/(6-a))*((4-2*a-d)/(6-a))^0.5)*(p^0.5))-K_H2O;
F(2) = (((2*d)/(6-a))/(((2-a-d)/(6-a))^(0.5)*((4-2*a-d)/(6-a))^(0.5)))-K_OH;
end

It returns

Undefined function or variable 'a'.

Error in eql_coef (line 2)
F(1) = (((2*a)/(6-a))/(((4-2*a-d)/(6-a))*((4-2*a-d)/(6-a))^0.5)*(p^0.5))-K_H2O;

I am pretty sure I have provided all the information I need and this is a syntac issue. Please help.

[Muhammad Usman]

your unknows are a & d and you are creating function for coefficients, that's seem wrong.

[Keith Contrera]

I thought you need to create a function to solve for variables a and d? How do you suggest I go about it?

Answer 1

[Muhammad Usman]

I think it will works like in this way:

function F = eql_coef(x)
% x(1) == a
% x(2) == d
K_H2O = 20.999;
K_OH = 1.0478;
p = 50;
F = [(((2*x(1))/(6-x(1)))/(((4-2*x(1)-x(2))/(6-x(1)))*((4-2*x(1)-x(2))/(6-x(1)))^0.5)*(p^0.5))-K_H2O;
    (((2*x(2))/(6-x(1)))/(((2-x(1)-x(2))/(6-x(1)))^(0.5)*((4-2*x(1)-x(2))/(6-x(1)))^(0.5)))-K_OH];
end

and then:

guess = [2;2];
options = optimoptions('fsolve','Display','iter');
[x,fval] = fsolve(@eql_coef,guess,options);

It produced results like:

                                         Norm of      First-order   Trust-region
 Iteration  Func-count     f(x)          step         optimality    radius
     0          3         850.247                           656               1
     1          6         712.211              1            127               1
     2          9         524.365              1           54.2               1
     3         12          506.84              1           39.3               1
     4         15         488.656              1           32.7               1
     5         16         488.656              1           32.7               1
     6         19         470.782           0.25           32.1            0.25
     7         20         470.782          0.625           32.1           0.625
     8         23         459.878        0.15625           31.6           0.156
     9         24         459.878       0.390625           31.6           0.391
    10         25         459.878      0.0976563           31.6          0.0977
    11         28         458.193      0.0244141           31.5          0.0244
    12         29         458.193      0.0610352           31.5           0.061
    13         30         458.193      0.0152588           31.5          0.0153
    14         33         457.931      0.0038147           31.5         0.00381
    15         34         457.931     0.00953674           31.5         0.00954
    16         37         457.766     0.00238419           31.5         0.00238
    17         38         457.766     0.00596046           31.5         0.00596
    18         41         457.664     0.00149012           31.5         0.00149
    19         42         457.664     0.00372529           31.5         0.00373
    20         43         457.664    0.000931323           31.5        0.000931
    21         46         457.648    0.000232831           31.5        0.000233
    22         47         457.648    0.000582077           31.5        0.000582
    23         50         453.253    0.000145519           31.4        0.000146
    24         51         453.253    0.000363798           31.4        0.000364
    25         54         453.247    9.09495e-05           31.4        9.09e-05
    26         55         453.247    0.000227374           31.4        0.000227
    27         58         453.243    5.68434e-05           31.4        5.68e-05
    28         59         453.243    0.000142109           31.4        0.000142
    29         62         453.241    3.55271e-05           31.4        3.55e-05
    30         65         453.235    8.88178e-05           31.4        8.88e-05
    31         66         453.235    0.000222045           31.4        0.000222
    32         67         453.235    5.55112e-05           31.4        5.55e-05
    33         68         453.235    1.38778e-05           31.4        1.39e-05
    34         69         453.235    3.46945e-06           31.4        3.47e-06
    35         72         453.234    8.67362e-07           31.4        8.67e-07
    36         73         453.234     2.1684e-06           31.4        2.17e-06
    37         74         453.234    5.42101e-07           31.4        5.42e-07
    38         75         453.234    1.35525e-07           31.4        1.36e-07
    39         76         453.234    3.38813e-08           31.4        3.39e-08
    40         79         453.234    8.47033e-09           31.4        8.47e-09
    41         80         453.234    2.11758e-08           31.4        2.12e-08
    42         83         453.234    5.29396e-09           31.4        5.29e-09
    43         84         453.234    1.32349e-08           31.4        1.32e-08
    44         87         453.234    3.30872e-09           31.4        3.31e-09
    45         88         453.234    8.27181e-09           31.4        8.27e-09
    46         91         453.234    2.06795e-09           31.4        2.07e-09
    47         92         453.234    5.16988e-09           31.4        5.17e-09
    48         95         453.234    1.29247e-09           31.4        1.29e-09
    49         96         453.234    3.23117e-09           31.4        3.23e-09
    50         99         453.234    8.07794e-10           31.4        8.08e-10
    51        100         453.234    2.01948e-09           31.4        2.02e-09
    52        103         453.234    5.04871e-10           31.4        5.05e-10
    53        104         453.234    1.26218e-09           31.4        1.26e-09
    54        107         453.234    3.15544e-10           31.4        3.16e-10
    55        108         453.234    7.88861e-10           31.4        7.89e-10
    56        109         453.234    1.97215e-10           31.4        1.97e-10
    57        112         453.234    4.93038e-11           31.4        4.93e-11
    58        113         453.234     1.2326e-10           31.4        1.23e-10
    59        116         453.234    3.08149e-11           31.4        3.08e-11
    60        117         453.234    7.70372e-11           31.4         7.7e-11
    61        120         453.234    1.92593e-11           31.4        1.93e-11
    62        121         453.234    4.81482e-11           31.4        4.81e-11
    63        124         453.234    1.20371e-11           31.4         1.2e-11
    64        125         453.234    3.00927e-11           31.4        3.01e-11
    65        128         453.234    7.52316e-12           31.4        7.52e-12
    66        129         453.234    1.88079e-11           31.4        1.88e-11
    67        132         453.234    4.70198e-12           31.4         4.7e-12

No solution found.

fsolve stopped because the relative size of the current step is less than the
default value of the step size tolerance squared, but the vector of function values
is not near zero as measured by the default value of the function tolerance.

<stopping criteria details>

Hope so it will helpful for you

Answer 2

[Bjorn Gustavsson]

You need to make the function eql_coef depend on the coefficients a and d, something like this:

function F = eql_coef(a_d,K_H2O, K_OH, p)
a = a_d(1);
d = a_d(2);
F(1) = (((2*a)/(6-a))/(((4-2*a-d)/(6-a))*((4-2*a-d)/(6-a))^0.5)*(p^0.5))-K_H2O;
F(2) = (((2*d)/(6-a))/(((2-a-d)/(6-a))^(0.5)*((4-2*a-d)/(6-a))^(0.5)))-K_OH;
end

Now you have a function that will depend on a and d, that fsolve can adjust. So you call fsolve like this:

a_d_guess = [2;2];
[a_d_solution] = fsolve(@(a_d) eql_coef(a_d,K_H2O, K_OH, p), a_d_guess)

HTH

[Keith Contrera]

This gives me imaginary solutions and I need real solutions. How would I fix this?

[Alex Sha]

Hi, I get the solution below:

a: 1.09358949887241
d: 0.427017129161841