Problem with fsolve: system of 6 nonlinear equations and more

1 回表示 (過去 30 日間)
Samuele Nigrelli
Samuele Nigrelli 2016 年 4 月 28 日
編集済み: Alex Sha 2019 年 12 月 2 日
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

サインインしてこの質問に回答する。

採用された回答

Walter Roberson
Walter Roberson 2016 年 4 月 28 日
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.
  1 件のコメント
Samuele Nigrelli
Samuele Nigrelli 2016 年 4 月 28 日
Exactly. I made a ridiculous mistake in calculating the avarages XP

サインインしてコメントする。

その他の回答 (3 件)

Alex Sha
Alex Sha 2019 年 12 月 2 日
編集済み: Alex Sha 2019 年 12 月 2 日
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 2016 年 4 月 28 日
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 2016 年 4 月 28 日
編集済み: John D'Errico 2016 年 4 月 28 日
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.

カテゴリ

Help Center および File Exchange で Equation Solving についてさらに検索

タグ

製品

Community Treasure Hunt

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

Start Hunting!

Translated by