Non linear system of equations

Hello,

I am trying to solve a nonlinear system with fsolve, but what i get are some negative numbers (they are concentrations so this has no sense) and i have two constants, which they are x(3) and x(4), but the fsolve changes their value. This is my function:

function f=ask3(x)
Hnh3=62;
Hco2=3.4*10^(-2);
Hhno3=2.1*10^5;
Hso2=1.2;
Hh2o2=7.5*10^(4);
Ho3=1.1*10^(-2);
Ka1=1.7*10^(-5);
Ks1=1.3*10^(-2);
Ks2=6.6*10^(-8);
Kc1=4.3*10^(-7);
Kc2=4.7*10^(-11);
K124=1.2*10^(-2);
K12=7.5*10^(7);
K22=13;
K0=2.4*10^(4);
K1=3.7*10^(5);
K2=1.5*10^(9);
Pnh3=5*10^(-9);
Pso2=20*10^(-9);
Phno3=3*10^(-9);
Pco2=330*10^(-6);
Po3=10^(-8);
Ph2o2=10^(-9);
Kw=10^(-14);
wL=10^(-7);
T=298;
R=0.082;
f(1)=x(5)+x(6)+2*x(7)+x(8)+2*x(9)+2*x(10)+x(11)+x(12)-x(1)-x(2)-3*x(3)-2*x(4);
f(2)=((x(22)*Pso2)/(1+x(22)*wL*R*T)-x(16));
f(3)=(((Hnh3*Pnh3*(1+((Ka1*x(1)/Kw))))))-x(13);
f(4)=(Hco2*Pco2*(1+(Kc1/x(1))+(Kc1*Kc2/(x(1)^2))))-x(14);
f(5)=x(10)+x(12)-x(17);
f(6)=Phno3*(3.2*10^(6)/x(1))-x(15);
f(7)=(Hso2*Pso2*(1+(Ks1/x(1))+((Ks1*Ks2)/(x(1)^2))))-x(22);
f(8)=(x(13)/(1+Kw/(ka1*x(1))))-x(2);
f(9)=x(5)*x(1)-Kw;
f(10)=(x(14)/(1+(x(1)/Kc1)+(Kc2/x(1))))-x(6);
f(11)=(x(14)/(1+(x(1)/Kc2)+((x(1)^2)/(Kc1*Kc2))))-x(7);
f(12)=(x(16)/(1+(x(1)/Ks1)+(Ks2/(x(1)))))-x(8);
f(13)=(x(16)/(1+(x(1)/ks2)+((x(1)^2)/(Ks1*Ks2))))-x(9);
f(14)=(K124*x(17)/(x(1)+K124))-x(10);
f(15)=(x(15)/(1+(x(1)/Kw)))-x(11);
f(16)=((K0*x(23)+K1*x(8)+K2*x(9))*x(18))-x(20);
f(17)=((K12*x(1)*x(19)*x(8))/(1+K22*x(1)))-x(21);
f(18)=((Ho3*Po3)/(1+Ho3*wL*R*T))-x(18);
f(19)=((Hh2o2*Ph2o2)/(1+Hh2o2*wL*R*T))-x(19);
f(20)=((x(1)*x(17))/(x(1)+K124))-x(12);
f(21)=((x(22)*Pso2)/(1+x(22)*wL*R*T))-x(23);
f(22)=x(3)-33.3*10^(-6);
f(23)=x(4)-2.5*10^(-6);
end

And this is what i get:

>> x0 = [10^(-5); 10^(-9); 33.3*10^(-6); 2.5*10^(-6); 10^(-9)*ones(20,1)];
>> options = optimoptions('fsolve','FunctionTolerance',1e-9,'StepTolerance',1e-9)
>> x=fsolve(@ask3,x0,options)
x =

  50.000000000000000
  50.000000000000000
  18.181823215108206
  18.181822907333242
  18.181823215108206
  18.181822907333242
  -4.545451270407963
  -4.545451314375789
 -31.818176784891797
 -31.818177092666755
  22.727274485516169
  22.727274221709031
  27.272725514483831
  27.272725778290965
   0.000000001000004
   0.000000001000000
   0.000000001000000
   0.000000001000000
   0.000000001000000
   0.000000001000000
   0.000000001000000
   0.000000001000000
   0.000000001000000
   0.000000001000000

If someone can tell what is going wrong, please help me!!!

Answer 1

[darova]

Hi. Take a look on your equations

f(8)=(x(13)/(1+Kw/(ka1*x(1))))-x(2);                        % must be Ka1. You have no ka1
f(13)=(x(16)/(1+(x(1)/ks2)+((x(1)^2)/(Ks1*Ks2))))-x(9);     % must be Ks2. You have no ks2

Also you have 24 initial conditions for 23 equations

[thanos tria]

I tried it but still no change, the main problem i think is why it changes my constants

[Torsten]

If they are constants, treat them as constants and not as solution variables.

That means that you give them fixed values in your code.

[darova]

Changed number of initial conditions, ka1 -> Ka1, ks2 -> Ks2:

x0 = [10^(-5); 10^(-9); 33.3*10^(-6); 2.5*10^(-6); 10^(-9)*ones(19,1)];
x = fsolve(@ask3,x0)

0,000869296201419250
0,458119098147945
3,33000000000000e-05        % x(3)
2,50000000000000e-06        % x(4)
1,15035588372175e-11
5,55000699660977e-09
3,00070710553115e-16
7,17822071442425e-15
5,44995556611017e-19
0,221522934115655
1,27038591283277e-10
0,0160474204294980
0,458119408147945
1,12255500072967e-05
11,0434164837404
7,65876570998069e-15
0,237570354545153
1,09999997043244e-10
6,33836740557945e-05
4,02296980578652e-19
2,93320819813032e-14
3,82938285499393e-07
7,65876570998069e-15

Answer 2

[Alex Sha]

multi-solutions:

1:

x1 0.00357502545408908

x2 1.88403841412371

x3 3.32997281546349E-5

x4 2.49981876993596E-6

x5 2.40013395038012E-10

x6 1.29843805677586E-8

x7 -1.10785702458008E-14

x8 -9.50278762729169E-14

x9 -1.54751427023488E-13

x10 0.821490444430723

x11 1.13507635946112E-11

x12 0.244737436865726

x13 1.88403872421433

x14 1.12213509281634E-5

x15 2.68529556593353

x16 -7.17800864234119E-14

x17 1.06622788143864

x18 1.10513400397213E-10

x19 6.33836740557957E-5

x20 -3.4405900132488E-14

x21 -1.54590979340091E-12

x22 1.11273723232722E-7

x23 -1.79622935418478E-14

2:

x1: 0.00668754292028887

x2: 3.52433511899043

x3: 3.32999972852601E-5

x4: 2.49999819166035E-6

x5: -5.09612458467555E-11

x6: 7.1628391318175E-10

x7: 3.31259628598224E-12

x8: -2.6856715490869E-12

x9: 5.36361000017502E-12

x10: 1.38080558762419

x11: 1.99020393453699E-12

x12: 0.769516385974346

x13: 3.52433542899133

x14: 1.12207214292351E-5

x15: 1.43550480565361

x16: -1.12501631394056E-12

x17: 2.15032197359763

x18: 1.10020823609164E-10

x19: 6.33836740558585E-5

x20: 1.12427655851418E-12

x21: -7.79878664942969E-11

x22: 7.06541871297799E-8

x23: -1.08615672436958E-12

[madhan ravi]

Post the code that you tried in MATLAB.