Problem with fsolve: system of 6 nonlinear equations and more

2 visualizzazioni (ultimi 30 giorni)
Samuele Nigrelli
Samuele Nigrelli il 28 Apr 2016
Modificato: Alex Sha il 2 Dic 2019
I'm trying to solve the following system of nonlinear equations with fsolve:
function F=equations(x)
F(1)=h_hot*(Tav_hot-Thot_membrane)-((Thot_in-Thot_out)*Mav_hot*CP_w)/Am;
F(2)=h_hot*(Tav_hot-Thot_membrane)-(lambda*Nd_l+hmembrane*(Thot_membrane-Tcold_membrane));
F(3)=(lambda*Nd_l+hmembrane*(Thot_membrane-Tcold_membrane))-Kw/dhelta3*(Tcold_membrane-Thot_coolingplate);
F(4)=Kw/dhelta3*(Tcold_membrane-Thot_coolingplate)-Kp/dhelta4*(Thot_coolingplate-Tcold_coolingplate);
F(5)=Kp/dhelta4*(Thot_coolingplate-Tcold_coolingplate)-h_cold*(Tcold_coolingplate-Tav_cold);
F(6)=h_cold*(Tcold_coolingplate-Tav_cold)-CP_w*Mcold_in*(Tcold_out-Tcold_in)/Am;
x0=[310 300 315 312 311 295] %[310 300 315 312 311 295];
options=optimset('Display','iter')
[x,fval,exit,output]=fsolve(@equations,x0,options)
equations_guess=equations(x0)
And it returns:
Norm of First-order Trust-region
Iteration Func-count f(x) step optimality radius
0 7 1.81204e+12 1.19e+10 1
1 14 1.7734e+12 1 1.55e+10 1
2 21 1.72466e+12 2.5 1.67e+10 2.5
3 22 1.72466e+12 2.5 1.67e+10 2.5
4 29 1.68914e+12 0.625 3.03e+10 0.625
5 30 1.68914e+12 1.5625 3.03e+10 1.56
6 37 1.64038e+12 0.390625 1.27e+11 0.391
7 38 1.64038e+12 0.976562 1.27e+11 0.977
8 39 1.64038e+12 0.244141 1.27e+11 0.244
9 46 1.60886e+12 0.0610352 5.91e+11 0.061
10 47 1.60886e+12 0.152588 5.91e+11 0.153
11 48 1.60886e+12 0.038147 5.91e+11 0.0381
12 55 1.5969e+12 0.00953674 1.43e+12 0.00954
13 56 1.5969e+12 0.00953674 1.43e+12 0.00954
14 63 1.57561e+12 0.00238419 1.59e+13 0.00238
15 64 1.57561e+12 0.00596046 1.59e+13 0.00596
16 65 1.57561e+12 0.00149012 1.59e+13 0.00149
17 66 1.57561e+12 0.000372529 1.59e+13 0.000373
18 73 1.56808e+12 9.31323e-05 7.44e+13 9.31e-05
19 74 1.56808e+12 0.000232831 7.44e+13 0.000233
20 75 1.56808e+12 5.82077e-05 7.44e+13 5.82e-05
21 82 1.56316e+12 1.45519e-05 1.02e+15 1.46e-05
22 83 1.56316e+12 3.63798e-05 1.02e+15 3.64e-05
23 84 1.56316e+12 9.09495e-06 1.02e+15 9.09e-06
24 91 1.56124e+12 2.27374e-06 7.19e+14 2.27e-06
25 92 1.56124e+12 2.27374e-06 7.19e+14 2.27e-06
26 99 1.56106e+12 5.68434e-07 6.46e+14 5.68e-07
27 106 1.56073e+12 5.68434e-07 7.09e+14 5.68e-07
28 113 1.56057e+12 5.68434e-07 6.43e+14 5.68e-07
29 120 1.56023e+12 5.68434e-07 7.04e+14 5.68e-07
30 127 1.5601e+12 5.68434e-07 6.43e+14 5.68e-07
31 134 1.55975e+12 5.68434e-07 7.01e+14 5.68e-07
32 141 1.55965e+12 5.68434e-07 6.44e+14 5.68e-07
33 148 1.55927e+12 5.68434e-07 6.99e+14 5.68e-07
34 155 1.5592e+12 5.68434e-07 6.46e+14 5.68e-07
35 162 1.55881e+12 5.68434e-07 6.96e+14 5.68e-07
36 163 1.55881e+12 5.68434e-07 6.96e+14 5.68e-07
37 170 1.55868e+12 1.42109e-07 7.35e+14 1.42e-07
38 177 1.55862e+12 1.42109e-07 7.28e+14 1.42e-07
39 184 1.55855e+12 1.42109e-07 7.36e+14 1.42e-07
40 191 1.55849e+12 1.42109e-07 7.3e+14 1.42e-07
41 198 1.55841e+12 1.42109e-07 7.38e+14 1.42e-07
42 205 1.55836e+12 1.42109e-07 7.32e+14 1.42e-07
43 212 1.55828e+12 1.42109e-07 7.4e+14 1.42e-07
44 219 1.55822e+12 1.42109e-07 7.34e+14 1.42e-07
45 226 1.55814e+12 1.42109e-07 7.42e+14 1.42e-07
46 233 1.55808e+12 1.42109e-07 7.36e+14 1.42e-07
47 240 1.558e+12 1.42109e-07 7.44e+14 1.42e-07
48 247 1.55794e+12 1.42109e-07 7.38e+14 1.42e-07
49 254 1.55785e+12 1.42109e-07 7.46e+14 1.42e-07
50 261 1.55779e+12 1.42109e-07 7.4e+14 1.42e-07
51 268 1.5577e+12 1.42109e-07 7.48e+14 1.42e-07
52 275 1.55764e+12 1.42109e-07 7.42e+14 1.42e-07
53 282 1.55755e+12 1.42109e-07 7.5e+14 1.42e-07
54 289 1.55749e+12 1.42109e-07 7.44e+14 1.42e-07
55 296 1.5574e+12 1.42109e-07 7.51e+14 1.42e-07
56 303 1.55733e+12 1.42109e-07 7.46e+14 1.42e-07
57 310 1.55724e+12 1.42109e-07 7.53e+14 1.42e-07
58 317 1.55717e+12 1.42109e-07 7.48e+14 1.42e-07
59 324 1.55707e+12 1.42109e-07 7.55e+14 1.42e-07
60 331 1.55701e+12 1.42109e-07 7.5e+14 1.42e-07
61 338 1.5569e+12 1.42109e-07 7.57e+14 1.42e-07
62 345 1.55683e+12 1.42109e-07 7.52e+14 1.42e-07
63 352 1.55673e+12 1.42109e-07 7.58e+14 1.42e-07
64 359 1.55666e+12 1.42109e-07 7.53e+14 1.42e-07
65 366 1.55655e+12 1.42109e-07 7.59e+14 1.42e-07
66 373 1.55648e+12 1.42109e-07 7.55e+14 1.42e-07
67 380 1.55637e+12 1.42109e-07 7.6e+14 1.42e-07
68 387 1.5563e+12 1.42109e-07 7.56e+14 1.42e-07
69 394 1.55618e+12 1.42109e-07 7.6e+14 1.42e-07
70 401 1.55612e+12 1.42109e-07 7.56e+14 1.42e-07
71 408 1.556e+12 1.42109e-07 7.6e+14 1.42e-07
72 415 1.55594e+12 1.42109e-07 7.56e+14 1.42e-07
73 422 1.55583e+12 1.42109e-07 7.59e+14 1.42e-07
74 429 1.55577e+12 1.42109e-07 7.55e+14 1.42e-07
75 436 1.55566e+12 1.42109e-07 7.57e+14 1.42e-07
76 443 1.55561e+12 1.42109e-07 7.54e+14 1.42e-07
77 450 1.55552e+12 1.42109e-07 7.55e+14 1.42e-07
78 457 1.55548e+12 1.42109e-07 7.52e+14 1.42e-07
79 464 1.55541e+12 1.42109e-07 7.52e+14 1.42e-07
80 471 1.55538e+12 1.42109e-07 7.49e+14 1.42e-07
81 478 1.55532e+12 1.42109e-07 7.49e+14 1.42e-07
82 485 1.5553e+12 1.42109e-07 7.47e+14 1.42e-07
83 492 1.55526e+12 1.42109e-07 7.46e+14 1.42e-07
84 499 1.55525e+12 1.42109e-07 7.44e+14 1.42e-07
85 506 1.55522e+12 1.42109e-07 7.44e+14 1.42e-07
86 507 1.55522e+12 1.42109e-07 7.44e+14 1.42e-07
87 514 1.55469e+12 3.55271e-08 7.77e+14 3.55e-08
88 515 1.55469e+12 8.88178e-08 7.77e+14 8.88e-08
89 522 1.55435e+12 2.22045e-08 8.05e+14 2.22e-08
90 523 1.55435e+12 5.55112e-08 8.05e+14 5.55e-08
91 530 1.55421e+12 1.38778e-08 8.28e+14 1.39e-08
92 531 1.55421e+12 3.46945e-08 8.28e+14 3.47e-08
93 538 1.5542e+12 8.67362e-09 8.24e+14 8.67e-09
94 545 1.55416e+12 8.67362e-09 8.28e+14 8.67e-09
95 546 1.55416e+12 2.1684e-08 8.28e+14 2.17e-08
96 553 1.55415e+12 5.42101e-09 8.3e+14 5.42e-09
97 554 1.55415e+12 1.35525e-08 8.3e+14 1.36e-08
98 561 1.55413e+12 3.38813e-09 8.32e+14 3.39e-09
99 562 1.55413e+12 8.47033e-09 8.32e+14 8.47e-09
100 569 1.55413e+12 2.11758e-09 8.32e+14 2.12e-09
101 576 1.55412e+12 5.29396e-09 8.27e+14 5.29e-09
102 583 1.55409e+12 5.29396e-09 8.33e+14 5.29e-09
103 584 1.55409e+12 1.32349e-08 8.33e+14 1.32e-08
104 591 1.55408e+12 3.30872e-09 8.3e+14 3.31e-09
105 598 1.55407e+12 8.27181e-09 8.23e+14 8.27e-09
106 605 1.55403e+12 8.27181e-09 8.31e+14 8.27e-09
Solver stopped prematurely.
fsolve stopped because it exceeded the function evaluation limit,
options.MaxFunEvals = 600 (the default value).
x =
1.0e+02 *
3.0994 - 0.0000i 3.0095 + 0.0000i 3.1319 + 0.0000i 3.1319 + 0.0000i 3.1148 - 0.0000i 2.9335 - 0.0000i
fval =
1.0e+05 *
-0.4949 + 0.1306i -0.4949 + 0.1306i -0.0036 - 0.0000i -0.5144 + 0.0000i -8.6309 - 0.0000i 8.9511 - 0.0000i
exit =
0
output =
iterations: 106
funcCount: 605
algorithm: 'trust-region-dogleg'
firstorderopt: 8.3052e+14
message: 'Solver stopped prematurely.
fsolve stopped because it exceeded the function evaluation limit,
options.MaxFunEvals = 60...'
equations_guess=
1.0e+05 *
-3.2452 -3.4664 0.2218 -0.4550 -8.7591 9.0376
Could anyone tell me what I have to do to get through the problem? Moreover I don't know why the solution x and fval become complex numbers. It's clear that increasing the number of iterations won't help!!
Thank you!!!!
The complete code is copied below, if anyone finds it necessary:
function solve()
clc
%Parameters
%modulo
%L=1; %lunghezza m
%W=0.5; %larghezza m
%spacer
%df= %filament diameter m
%deltha1=3; %thickness m
%es= %porosity
%THspacer= %angle porosity
%hot channel
deltha1=0.003; %thickness m
%dh= %equivalent diameter m
Apax=deltha1*0.03; %trasverse surface m^2
Deq=2*deltha1; %equivalent diameter m
%membrane
epsilon=0.8; %porosity
Am=0.0417; %surface m^2
dp=2.00E-07; %pore diameter m
dhelta2 =2.50E-04; %thickness m
Km= 0.25; %conductivity W/mK
tau= 1.5; %tortuosity
K1= epsilon/tau;
K0= (2/3)*K1*(dp/2);
%airgap
dhelta3= 0.003; %thickness m
%Kag= %conductivity water-air W/mK
%coolant plate
dhelta4= 7.00E-05; %thickness m
Kp=0.2; %conductivity W/mK
%Antoine law log(10)P[mmHg]=A-B/(C+T[°C])
A=8.07131; B=1730.63; C=233.426;
%other parameters
lambda=2200000; %latent heat of vapor j/kg
Patm= 101325; %Pa
CP_w= 4186; %Cp water j/kgK
CP_vap= 1900; %Cp vap j/kgK
MU_cold= 8.55E-04; %viscosity cold Pa*s
MU_hot= 3.65E-04; %viscosity hot Pa*s
PM= 0.018; % KG/mol water
R= 8.314; %universal constant Pa*m^3/(mol*k)
Kw= 0.63; %conductivity water W/mK
Kair= 0.03; %conductivity air W/mK
D12= 2.20E-05; %diffusivity water-air m^2/s
PR_hot= CP_w*MU_hot/Kw; %prandtl hot channel
PR_cold= CP_w*MU_cold/Kw; %prandtl cold channel
%boundary condition
Mhot_in= 0.02 ; %Feed mass in kg/s
Thot_in= 323.2 ; %Feed temperature in K
Mcold_in= 0.02 ; %Coolant mass in kg/s
Tcold_in= 291.1 ; %Coolant temperature in K
%v_hot=Mhot_in/(Apax*1000); %speed hot feed
%resolution
x0=[310 300 315 312 311 295]; %[310 300 315 312 311 295]
options=optimset('Display','iter')
[x,fval,exit,output]=fsolve(@equations,x0,options)
equations_guess=equations(x0)
function F=equations(x)
%J=x(1); %heat flux W/m^2
Thot_out= x(1); %Feed temperature out K
Tcold_out= x(2); %Coolant temperature out K
Thot_membrane= x(3); %temperature on the membrane hot side K
Tcold_membrane= x(4); %temperature on the membrane permeate side K
Thot_coolingplate= x(5); %temperature on cooling plate hot K
Tcold_coolingplate= x(6); %temperature on cooling plate cold K
Tav_hot= (Thot_in-Thot_out)/2; %temperature avarage in-out hot K
Tav_cold= (Tcold_out-Tcold_in)/2; %temperature avarage in-out cold K
Tav_m= (Thot_membrane-Tcold_membrane)/2; %temperature avarage membrane K
Psat_hot=133.3*10^(A-(B/(C+Thot_membrane-273))); %saturation pressure
Psat_cold=133.3*10^(A-(B/(C+Tcold_membrane-273))); %saturation pressure
hmembrane=1/dhelta2*(epsilon*Kair+(1-epsilon)*Km); %membrane transfer heat coefficent
Pair_Hot_membrane= Patm-Psat_hot;
Pair_Cold_membrane= Patm-Psat_cold;
Pair_ml=(Pair_Hot_membrane-Pair_Cold_membrane)/log(Pair_Hot_membrane/Pair_Cold_membrane);
Dknudsen=K0*(8*R*Tav_m/(pi*PM))^(1/2); %knudsen diffusivity
Dstagnant=D12*Patm/Pair_ml; %stagnant diffusivity
Koverall=1/(R*Tav_m)*(1/((1/Dknudsen+1/Dstagnant))); %membrane transfer mass coefficent
Nd=Koverall/dhelta2*(Psat_hot-Psat_cold); %distillate flux mol/m^2*s
Nd_l=Nd*PM; %distillate flux l/m^2*s o kg/m^2*s
Mhot_out= Mhot_in-Nd_l*Am; %mass hot out kg/s
Mav_hot= (Mhot_in-Mhot_out)/2; %mass avarage hot kg/s
Re_hot=Mav_hot*Deq/(Apax*MU_hot); %Reinolds hot channel
Re_cold=Mcold_in*Deq/(Apax*MU_cold); %Reinolds cold channel
h_hot=((8.4577*(Re_hot^0.1652)+0.7214*(Re_hot^0.5155))/2)*Kw/Deq; %transfer heat coeffcient hot W/m^2K
h_cold=((8.4577*(Re_cold^0.1652)+0.7214*(Re_cold^0.5155))/2)*Kw/Deq; %transfer heat coeffcient cold W/m^2K
F(1)=h_hot*(Tav_hot-Thot_membrane)-((Thot_in-Thot_out)*Mav_hot*CP_w)/Am;
F(2)=h_hot*(Tav_hot-Thot_membrane)-(lambda*Nd_l+hmembrane*(Thot_membrane-Tcold_membrane));
F(3)=(lambda*Nd_l+hmembrane*(Thot_membrane-Tcold_membrane))-Kw/dhelta3*(Tcold_membrane-Thot_coolingplate);
F(4)=Kw/dhelta3*(Tcold_membrane-Thot_coolingplate)-Kp/dhelta4*(Thot_coolingplate-Tcold_coolingplate);
F(5)=Kp/dhelta4*(Thot_coolingplate-Tcold_coolingplate)-h_cold*(Tcold_coolingplate-Tav_cold);
F(6)=h_cold*(Tcold_coolingplate-Tav_cold)-CP_w*Mcold_in*(Tcold_out-Tcold_in)/Am;
end
save
end

Accedi per rispondere a questa domanda.

Risposta accettata

Walter Roberson
Walter Roberson il 28 Apr 2016
Your Thot_membrane is becoming less than your Tcold_membrane which gives a negative Tav_m which leads to
Dknudsen=K0*(8*R*Tav_m/(pi*PM))^(1/2)
taking the square root of a negative number.

Più risposte (3)

Alex Sha
Alex Sha il 2 Dic 2019
Modificato: Alex Sha il 2 Dic 2019
a stable results:
thot_membrane 470.276797155953
tcold_membrane 470.271756744994
thot_out -727.941719280904
thot_coolingplate 75.7496078549212
tcold_coolingplate 46.7522299114363
tcold_out 332.366488985555
Torsten
Torsten il 28 Apr 2016
My guess is that Re_hot and/or Re_cold become negative in the course of the iteration process.
Best wishes
Torsten.
John D'Errico
John D'Errico il 28 Apr 2016
Modificato: John D'Errico il 28 Apr 2016
It just might be there are no solutions. Not every problem has a real solution. I set up your problem as a symbolic one.
res = vpasolve(F{:})
res =
Tcold_coolingplate: [1x1 sym]
Tcold_membrane: [1x1 sym]
Tcold_out: [1x1 sym]
Thot_coolingplate: [1x1 sym]
Thot_membrane: [1x1 sym]
Thot_out: [1x1 sym]
res.Tcold_out
ans =
302.43921796117817410514164884414 - 121.66445227286058883652318399616i
res.Tcold_coolingplate
ans =
12.846591463673943345342329712164 - 137.83785789745060515861846007314i
res.Tcold_membrane
ans =
129.22141141302985703547817144583 - 1386.4847025004829566433917220377i
res.Tcold_out
ans =
302.43921796117817410514164884414 - 121.66445227286058883652318399616i
res.Thot_coolingplate
ans =
20.814499490015498684163926700915 - 223.32974711805878199525640581548i
res.Thot_membrane
ans =
137.91215125253110323480833978171 + 83.998157315424176761677478609752i
res.Thot_out
ans =
- 55.34925693520510435933273846792 - 23.965610390231809575830847603116i
All it could find are complex solutions. Lots of fractional exponents of variables in there, so I'm not at all surprised at seeing complex results in one spot or another. And once one variable has an imaginary part, it will propagate.

Categorie

Scopri di più su Spline Postprocessing in Help Center e File Exchange

Tag

Prodotti

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by