应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 程序语言 » MATLAB/Mathematica » matlab实现牛顿迭代法求非线性方程组程序运行出现问题
91/1返回列表
查看: 1480  |  回复: 8
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

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;

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

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:
1楼 2015-04-20 14:44:55
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

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的回帖

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的回帖

dingd

铁杆木虫 (职业作家)

【答案】应助回帖

感谢参与，应助指数 +1
Matlab感觉好复杂的样子。1stOpt求解看看：
CODE:
Constant T=318.2, rhomcA=0.34, CA=3.26, EAA=64, EAs=-1385, rhomcB=0.98, CB=4.53, EBB=82, EBs=-1690;
LoopConstant 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];
Function
ln(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;
ln(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;

Loop constant rhob        Loop constant ya        Loop constant yb        x1        x2
0.0047        0.89        0.11        0.349788429554512        0.0364739534111791
0.009        0.886        0.114        0.460684483692858        0.0489828447137933
0.0177        0.884        0.116        0.580741252753925        0.0617304175836032
0.0264        0.883        0.117        0.649820821824047        0.0690991828120868
0.0367        0.877        0.123        0.699996476771262        0.0783646455649231
0.0466        0.871        0.129        0.732280709438354        0.0863910296157365
0.0566        0.866        0.134        0.756187405093584        0.0931836288289366
0.0674        0.861        0.139        0.775229538588392        0.0997840413222308
0.0787        0.856        0.144        0.789855870165985        0.106177522100049
0.089        0.852        0.148        0.800097262910873        0.111365586073744
赞一下(1人) 回复此楼
6楼2015-04-20 17:32:40
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

cooooldog

铁杆木虫 (著名写手)

ส็็็

【答案】应助回帖

★ ★ ★ ★ ★ ★ ★ ★ ★ ★
感谢参与，应助指数 +1
xwl196: 金币+10, ★★★很有帮助 2015-04-20 21:03:13
1. 建议使用Matlab内部函数fsolve来做,不要自己写了;

2. 这类问题用符号计算求Jacobian意义不大,通常更麻烦而已,差分法效果通常更好;

3. 好好研究下fsolve里面的选项和用法
赞一下(1人) 回复此楼
ส็็็็็็็็็็็็็็็็็็็็
7楼2015-04-20 19:55:59
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

xwl196

银虫 (小有名气)
引用回帖:
6楼: Originally posted by dingd at 2015-04-20 17:32:40
Matlab感觉好复杂的样子。1stOpt求解看看：

Constant T=318.2, rhomcA=0.34, CA=3.26, EAA=64, EAs=-1385, rhomcB=0.98, CB=4.53, EBB=82, EBs=-1690;
LoopConstant rhob=,
             yA=,
             y ...

大神用的是1stOpt那个版本？
赞一下 回复此楼
8楼2015-04-20 21:02:48
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

xwl196

银虫 (小有名气)
引用回帖:
6楼: Originally posted by dingd at 2015-04-20 17:32:40
Matlab感觉好复杂的样子。1stOpt求解看看：

Constant T=318.2, rhomcA=0.34, CA=3.26, EAA=64, EAs=-1385, rhomcB=0.98, CB=4.53, EBB=82, EBs=-1690;
LoopConstant rhob=,
             yA=,
             y ...

1stopt1.5版本的计算不了，高版本的没有版权，悲催了
赞一下 回复此楼
9楼2015-04-20 21:11:40
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 xwl196 的主题更新
91/1返回列表
普通表情 龙 兔 虎 猫
信息提示
关闭
请填处理意见
关闭