24小时热门版块排行榜    

查看: 1536  |  回复: 8
当前只显示满足指定条件的回帖,点击这里查看本话题的所有回帖

xwl196

银虫 (小有名气)

[求助] matlab实现牛顿迭代法求非线性方程组程序运行出现问题 已有3人参与

已知非线性方程组如下:
log(rhomcA*x1*(1-yA*rhob/rhomcA-yB*rhob/rhomcB)/(yA*rhob*(1-x1-x2)))+EAA/T*(7*x1-8*yA*rhob/rhomcA)+(EAA*EBB)^0.5/T*(7*x2-8*yB*rhob/rhomcB)+EAs/T=0
f2=log(rhomcB*x2*(1-yA*rhob/rhomcA-yB*rhob/rhomcB)/(yB*rhob*(1-x1-x2)))+EBB/T*(7*x2-8*yB*rhob/rhomcB)+(EAA*EBB)^0.5/T*(7*x1-8*yA*rhob/rhomcA)+EBs/T=0
首先建立函数fun
function f=fun(x);
%定义非线性方程组如下
%变量x1 x2
%函数f1 f2
T=318.2;     %(K)
rhomcA=0.34; %(g/cm3)
CA=3.26;     %(mmol/g)
EAA=64;      %(K)
EAs=-1385;   %(K)
rhomcB=0.98; %(g/cm3)
CB=4.53;     %(mmol/g)
EBB=82;      %(K)
EBs=-1690;   %(K)
rhob=[0.0047 0.0090 0.0177 0.0264 0.0367 0.0466 0.0566 0.0674 0.0787 0.0890]';
yA=[0.89 0.886 0.884 0.883 0.877 0.871 0.866 0.861 0.856 0.852]';
yB=[0.11 0.114 0.116 0.117 0.123 0.129 0.134 0.139 0.144 0.148]';
syms x1 x2 x3 x4
f1=log(rhomcA*x1.*(1-yA.*rhob/rhomcA-yB.*rhob/rhomcB)./(yA.*rhob.*(1-x1-x2)))+EAA/T*(7*x1-8*yA.*rhob/rhomcA)+(EAA*EBB)^0.5/T*(7*x2-8*yB.*rhob/rhomcB)+EAs/T;
f2=log(rhomcB*x2.*(1-yA.*rhob/rhomcA-yB.*rhob/rhomcB)./(yB.*rhob.*(1-x1-x2)))+EBB/T*(7*x2-8*yB.*rhob/rhomcB)+(EAA*EBB)^0.5/T*(7*x1-8*yA.*rhob/rhomcA)+EBs/T;
%f3=nexA-2*CA*(x1-x3*rhob/rhomcA);
%f4=nexB-2*CB*(x2-x4*rhob/rhomcB);
f=[f1 f2];

建立函数dfun
function df=dfun(x);
%用来求解方程组的雅克比矩阵储存在dfun中
f=fun(x);
df=[diff(f,'x1');diff(f,'x2')];
df=conj(df');

编程牛顿法求解非线性方程组将newton.m保存到工作路径中:
function x=newton(x0,eps,N);
con=0;
%其中x0为迭代初值,eps为精度要求N为最大迭代步数,con用来记录结果是否收敛
for i=1:N;
    f=subs(fun(x0),{'x1' 'x2'},{x0(1) x0(2)});
    df=subs(dfun(x0),{'x1' 'x2'},{x0(1) x0(2)});
    x=x0-f./df;
    for j=1:length(x0);
        il(i,j)=x(j);
    end
    if norm(x-x0)<eps
       con=1;
       break;
    end
    x0=x;
end

%以下是将迭代过程写入txt文档文件名为iteration.txt

fid=fopen('iteration.txt','w');
fprintf(fid,'iteration');
for j=1:length(x0)
  fprintf(fid,'         x%d',j);
end
for j=1:i
fprintf(fid,'\n%6d     ',j);
for k=1:length(x0)
fprintf(fid,' %10.6f',il(j,k));
end
end
if con==1
  fprintf(fid,'\n计算结果收敛!');
end
if con==0
fprintf(fid,'\n迭代步数过多可能不收敛!');
end
fclose(fid);


在matlab中命令窗口中输入:
newton([0.8 0.8],0.00001,20)

结果出现:
??? Error using ==> rdivide
Matrix dimensions must agree.

Error in ==> newton at 7
    x=x0-f./df;

检查了几遍程序没发现问题出现在哪,请高手指点
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐,对您同样有帮助:

已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

xwl196

银虫 (小有名气)

引用回帖:
4楼: Originally posted by mylifeljy at 2015-04-20 15:42:05
哈哈,这个也没什么成文或不成文的规定。
综合考虑你的财富与问题,可以发现你是个勤俭节约的好妹子。。。玩笑而已,不必当真~
...

按着你说的可能出现的问题,我对程序又改了改,发现还是没结果,然后捋了捋思路,问题可能出现在如果我按着上述程序计算,即使能得出结果,对应的x1和x2也应该是数,而不是数组。
我把我要解决的问题给你描述一下:
以下参数已知:
T=318.2;     
rhomcA=0.34;
CA=3.26;     
EAA=64;      
EAs=-1385;   
rhomcB=0.98;
CB=4.53;     
EBB=82;      
EBs=-1690;   
rhob=[0.0047 0.0090 0.0177 0.0264 0.0367 0.0466 0.0566 0.0674 0.0787 0.0890]';
yA=[0.89 0.886 0.884 0.883 0.877 0.871 0.866 0.861 0.856 0.852]';
yB=[0.11 0.114 0.116 0.117 0.123 0.129 0.134 0.139 0.144 0.148]';
已知公式如下:
log(rhomcA*x1*(1-yA*rhob/rhomcA-yB*rhob/rhomcB)/(yA*rhob*(1-x1-x2)))+EAA/T*(7*x1-8*yA*rhob/rhomcA)+(EAA*EBB)^0.5/T*(7*x2-8*yB*rhob/rhomcB)+EAs/T=0
log(rhomcB*x2*(1-yA*rhob/rhomcA-yB*rhob/rhomcB)/(yB*rhob*(1-x1-x2)))+EBB/T*(7*x2-8*yB*rhob/rhomcB)+(EAA*EBB)^0.5/T*(7*x1-8*yA*rhob/rhomcA)+EBs/T=0
想使用newton迭代法求x1和x2,要求的x1和x2应该分别为10×1数据
请大神帮帮忙呗,这个问题困扰我一周了,matlab初学者,一遇到循环就蒙
5楼2015-04-20 16:33:16
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 9 个回答

mylifeljy

禁虫 (正式写手)

感谢参与,应助指数 +1
本帖内容被屏蔽

2楼2015-04-20 15:18:07
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

xwl196

银虫 (小有名气)

引用回帖:
2楼: Originally posted by mylifeljy at 2015-04-20 15:18:07
妹子,你也太节俭了吧,都这么有钱了才出10个金币来解决这么冗长的问题!
如matlab报错提示,你的第一个问题是出在newton函数中, x=x0-f./df这句, 因为subs计算后f为一个10X2大小的矩阵,而df为一个2X20大小的矩阵 ...

10个金币少啊,不懂规矩,莫怪
多谢!
3楼2015-04-20 15:33:42
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

mylifeljy

禁虫 (正式写手)

本帖内容被屏蔽

4楼2015-04-20 15:42:05
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
信息提示
请填处理意见