非线性方程组求解问题（关于TolFun设置问题）

在使用fsolve求解非线性方程时，

得到的提示是

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the default value of the function tolerance, and
the problem appears regular as measured by the gradient.

<stopping criteria details>


Equation solved. The sum of squared function values, r = 1.810346e-19, is less than
sqrt(options.TolFun) = 1.000000e-03. The relative norm of the gradient of r, 3.334577e-10,
is less than options.TolFun = 1.000000e-06.


于是我修改了TolFun的范围至e-20

提示变成了
Equation solved, inaccuracy possible.

The vector of function values is near zero, as measured by the selected value
of the function tolerance. However, the last step was ineffective.

<stopping criteria details>


fsolve stopped because the sum of squared function values, r, changed by 3.037304e-25
relative to its initial value; this is less than max(options.TolFun^2,eps) = 2.220446e-16.
r = 9.693523e-27, is less than sqrt(options.TolFun) = 1.000000e-10.

Optimization Metric                                                Options
relative change r =   3.04e-25                 max(TolFun^2,eps) =   2e-16 (selected)
r =   9.69e-27                                     sqrt(TolFun) =  1.0e-10 (selected)

看提示是说，TolFun^2，eps的范围也偏大了，
那么再哪里调整这个值呢？
有可能我说的存在问题，
谢谢大家了，麻烦帮我看看应该怎么修改。

以下是代码

主代码部分：

1. function main
   
2. clear all
   
3. clc
   
4. global n G dl M M1 b1 b2 C4 C6 C7 C8 C9 a
   
5. T=273.2+112.13;                                 %反应器内平均温度,K
   
6. RMI=1.1;                                        %氯化氢乙炔分子比
   
7. u=0.062*3.5;                                    %气体线速度,m/s
   
8. C_A=exp(1.1835+1.055e-3*T)*4.184e3;             %乙炔热容,J/kg/K
   
9. C_H=0.192*4.184e3;                              %氯化氢热容,J/Kg/K
   
10. C_V=exp(2.238+1.256e-3*T)*4.184*1000/62.5;      %氯乙烯热容,J/Kg/K
    
11. L_A=exp(9.6084+1.413*log(T))*0.1;               %乙炔导热系数,W/m
    
12. L_H=exp(-7.574+0.9769*log(T))*0.1;              %氯化氢导热系数,W/m/K
    
13. L_V=0.018;                                      %氯乙烯导热系数,W/m/K
    
14. r_A=0.0115*(T/353.2)^0.97*0.001;                %乙炔黏度,Pa.s
    
15. r_H=0.0165*(T/353.2)^0.95*0.001;                %氯化氢黏度,pa.s
    
16. r_V=8.601e-3*(T/250.1)^0.97*0.001;              %氯乙烯黏度,Pa.s
    
17. rou_A=0.9130;                                   %乙炔密度,Kg/m3
    
18. rou_H=1.2817;                                   %氯化氢密度,Kg/m3
    
19. rou_V=2.1946;                                   %氯乙烯密度,Kg/m3
    
20. x_A=0.5;
    
21. C_av=(C_A*(1.0-x_A)+C_H*(RMI-x_A)+C_V*x_A)/(RMI+1.0-x_A);           %平均热容,J/Kg/K
    
22. L_av=(L_A*(1.0-x_A)+L_H*(RMI-x_A)+L_V*x_A)/(RMI+1.0-x_A);           %平均导热系数,W/m/K
    
23. r_av=(r_A*(1.0-x_A)+r_H*(RMI-x_A)+r_V*x_A)/(RMI+1.0-x_A);           %平均黏度,Pa.s
    
24. rou_av=(rou_A*(1.0-x_A)+rou_H*(RMI-x_A)+rou_V*x_A)/(RMI+1.0-x_A);   %平均密度,Kg/m3
    
25. re=0.044*u*rou_av/r_av;                                             %雷诺数
    
26. pr=C_av*r_av/L_av;                                                  %普朗特数
    
27. l_er=L_av*(10.3+0.037*re*pr);                                       %床层有效导热系数,W/(mK)
    
28. dp=0.004327;                                                        %催化剂颗粒直径,m
    
29. hw=(2.58*re^(1/3)*pr^(1/3)+0.094*re^0.8*pr^0.4)*L_av/dp;            %床层与器壁间导热系数,W/(m2*K)
    
30. 
    
31. %网格和步长
    
32. n=5+1;                  %径向有6个点,包括圆心
    
33. m=268;                  %？？？m是什么
    
34. dr=0.0044;              %径向微元长度
    
35. dl=0.0112;              %轴向微元长度
    
36. 
    
37. %参数
    
38. Cp=1.206e3;             %反应混合物比热容,J/(kg*K)
    
39. dH=-109e3;              %反应热,J/mol
    
40. rou_b=610;              %催化剂堆积密度,kg/m^3
    
41. G=u*rou_av;             %总质量速率,kg/(m^2.s)
    
42. v=u*3600/3/2.1;         %乙炔空速,hr-1
    
43. C_A0=20.7;              %进口乙炔浓度,mol/m3
    
44. Er=dp*u/10;             %径向混合扩散系数,m2/s
    
45. Tw=371.2;               %管壁温度,K
    
46. 
    
47. %方程的系数
    
48. a1=l_er/(G*Cp);         %(3-1)方程中的a1,b1
    
49. b1=-dH*rou_b/G/Cp;      
    
50. a2=Er/u;                %(3-2)方程中的a1,b1
    
51. b2=rou_b/u/C_A0;
    
52. M=a1*dl/dr^2;
    
53. M1=a2*dl/dr^2;
    
54. C1=4*dr*hw/(-l_er);     %(3-25)(3-26)中的系数c1-c9 ？？？c1计算公式
    
55. C2=(1+C1*M*(1-1/2/n)/4/(1+M))/(1-C1*M*(1+1/2/n)/4/(1+M));
    
56. C3=C1/(1+M)/(1-C1*M*(1+1/2/n)/4/(1+M));
    
57. C4=b1*C1*dl/4/(1+M)/(1-C1*M*(1+1/2/n)/4/(1+M));
    
58. C5=C1*Tw/(1-C1*M*(1+1/2/n)/4/(1+M));
    
59. C6=M*(1-1/2/n+C2*(1+1/2/n))/2/(1+M);
    
60. C7=(M*C3*(1+1/2/n))/2/(1+M)+(1-M)/(1+M);
    
61. C8=(M*C4*(1+1/2/n)/2/(1+M)+b1*dl/2)/(1+M);
    
62. C9=M*C5*(1+1/2/n)/2/(1+M);
    
63. 
    
64. %初始变量,反应器入口
    
65. a=0.209;                %？？？a=0.209
    
66. T(1:n-1)=362.8;         %进口温度
    
67. T(n)=371.2;             %管壁温度
    
68. x(1:n)=0;               %入口氯乙烯浓度
    
69. XO=[T;x];
    
70. Tresult(1,:)=XO(1,:);
    
71. xresult(1,:)=XO(2,:);
    
72. 
    
73. %i：径向，j：轴向
    
74. %计算轴向起点处的F(i,j), G(i,j) 和 rc(i,j)
    
75. %计算对应于j=1时的F(i,1), G(i,1)
    
76. FGrc(T,x);
    
77. 
    
78. % 沿轴向移动
    
79.     options=optimset('tolx',1e-20);      
    
80.     options=optimset(options,'tolfun',1e-20);
    
81. for j=1:m-1
    
82.     j=j+1;
    
83.     X=fsolve(@TxEquations,XO,options);  %求解非线性方程组，代入初值为X0
    
84.     T=X(1,:);
    
85.     x=X(2,:);
    
86.     Tresult(j,:)=T;
    
87.     xresult(j,:)=x;
    
88.    
    
89.     %轴向平均温度，转化率用梯形面积法计算？？？没有中点温度，且壁面温度多了一份，是否认为T(1)=T(6)
    
90.     xa(j,:)=sum(2*xresult(j,2)+4*xresult(j,3)+6*xresult(j,4)+8*xresult(j,5)+5*xresult(j,6))/25;%径向均温
    
91.     Ta(j,:)=sum(2*Tresult(j,2)+4*Tresult(j,3)+6*Tresult(j,4)+8*Tresult(j,5)+5*Tresult(j,6))/25;%径向平均转化率
    
92.     %计算F(i,j)和G(i,j)以及反应速度re(i,j)供下次迭代使用
    
93.     FGrc(T,x)
    
94. end
    
95. 
    
96. l=dl.*[0:j-1];
    
97. r=dr.*[0:n-1];
    
98. l=l';
    
99. %结果作图
    
100. surf(r,l,Tresult)   %反应管轴径向温度分布
     
101. xlabel('r(m)')
     
102. ylabel('l(m)')
     
103. zlabel('T(K)')
     
104. figure
     
105. plot(l,Ta)          %径向均温轴向分布(模拟出的均温有很大波动)
     
106. xlabel('l(m)')
     
107. ylabel('T_a_v(K)')
     
108. figure
     
109. plot(l,xa)          %径向转化率轴向分布
     
110. xlabel('l(m)')
     
111. ylabel('j_a_v')
     
112. figure
     
113. surf(r,l,xresult)   %反应管轴径向转化率分布
     
114. xlabel('r(m)')
     
115. ylabel('l(m)')
     
116. zlabel('x')

复制代码

函数：FGrc

1. function f=FGrc(T,x)
   
2. % 该脚本功能：计算某特定轴向距离时的F(i),G(i)
   
3. % i：径向，j：轴向
   
4. % 以下用T1代表执行第一次
   
5. % T1：j=1，即为入口处。
   
6. global F G rc n dl M M1 b1 b2 C6 C7 C8
   
7. rc=ReactionRate(T,x);                                                           %计算反应速率
   
8. 
   
9. % T1：计算入口处的圆心点数据
   
10. F(1)=((1-2*M)*T(1)+2*M*T(2)+b1*dl/2*rc(1))/(1+2*M);                             %(3-20)
    
11. G(1)=((1-2*M1)*x(1)+2*M1*x(2)+b2*dl/2*rc(1))/(1+2*M1);                          %(3-21)
    
12. 
    
13. % T1：计算入口处的非壁面的其他点数据
    
14. i=(2:5);
    
15.     var1=1-0.5./(i-1);
    
16.     var2=1+0.5./(i-1);
    
17.     F(i)=(M/2*(var1.*T(i-1)+var2.*T(i+1))+(1-M)*T(i)+b1*dl/2*rc(i))/(M+1);      %(3-12)
    
18.     G(i)=(M1/2*(var1.*x(i-1)+var2.*x(i+1))+(1-M1)*x(i)+b2*dl/2*rc(i))/(M1+1);   %(3-13)
    
19.    
    
20. % T1：计算入口处的壁面点数据
    
21. F(n)=C6*T(n-1)+C7*T(n)+C8*rc(n);                                                %(3-28)
    
22. G(n)=(M1*x(n-1)+(1-M1)*x(n)+b2*dl/2*rc(n))/(1+M1);                              %(3-23)

复制代码

方程组

1. function f=TxEquations(X)
   
2. % 该脚本功能：求解方程组
   
3. % i：径向，j：轴向
   
4. global n F G rc dl M M1 b1 b2 C6 C8 C9
   
5. T=X(1,:);
   
6. x=X(2,:);
   
7. % 圆心点温度组成，用到了轴向上一点F,G。与轴向该点近圆心点连立方程。
   
8. fT(1)=F(1)+(2*M*T(2)+b1/2*dl*rc(1))/(1+2*M)-T(1);       %(3-18)
   
9. fx(1)=G(1)+(2*M1*x(2)+b2/2*dl*rc(1))/(1+2*M1)-x(1);     %(3-19)
   
10. 
    
11. for i=2:n-1
    
12.     var1(i)=(1-0.5/(i-1));
    
13.     var2(i)=(1+0.5/(i-1));
    
14.     fT(i)=F(i)+(M/2*(var1(i)*T(i-1)+var2(i)*T(i+1))+b1/2*dl*rc(i))/(M+1)-T(i);   %(3-10)
    
15.     fx(i)=G(i)+(M1/2*(var1(i)*x(i-1)+var2(i)*x(i+1))+b2/2*dl*rc(i))/(M1+1)-x(i); %(3-11)
    
16. end
    
17. fT(n)=F(n)+C6*T(n-1)+C8*rc(n)-C9-T(n);                  %(3-27)
    
18. fx(n)=G(n)+(M1*x(n-1)+b2/2*dl*rc(n))/(1+M1)-x(n);       %(3-22)
    
19. f=[fT;fx];

复制代码

函数：ReactionRate

1. function f=ReactionRate(T,x)
   
2. % 该脚本功能：计算不同温度，组成时的反应速度
   
3. % ？？？a是什么？？？
   
4. global a
   
5. P=1.48;                             %固定床入口处压力,atm
   
6. RMI=1.1;                            %HCL,C2H2摩尔比
   
7. C_0=0.9935/(1+RMI);                 %乙炔进口相对摩尔分率
   
8. P_AO=P*C_0;                         %乙炔分压,atm
   
9. R=1.987;                            %气体常数,cal/mol/K
   
10. yita=57./(T-300)-0.12*exp(0.1.*x);  %催化剂有效系数
    
11. k=25.6.*yita.*exp(-4065/R./T);      %反应速率常数
    
12. K_H=3.31e6*exp(-10630/R./T);        %HCL的吸附平衡常数
    
13. %ra的表达式
    
14. f=(P_AO^2*a*k.*K_H/3.6).*(1-x).*(RMI-x)./((P-P_AO.*x).^2+P_AO*K_H.*(P-P_AO.*x).*(RMI-x));
    

复制代码

你好，我遇到的问题和你的相似，能留个qq交流下么？