24小时热门版块排行榜    

查看: 1769  |  回复: 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的回帖

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
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里面的选项和用法
ส็็็็็็็็็็็็็็็็็็็็
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 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 欢迎发来filecode的Mz6后的代码验证其规律 +52 医学老男孩 2026-08-13 112/5600 2026-08-19 07:36 by wuhaigang
[基金申请] 重要消息,中午系统在维护 +10 yuleib84 2026-08-18 11/550 2026-08-19 07:30 by 羽毛枫f
[基金申请] 什么时候开奖? +4 CrisMessi 2026-08-18 4/200 2026-08-19 00:53 by gatelove
[基金申请] 咱们一起用铁证分析2026国家社科基金中标与否 +7 启萌科技 2026-08-12 27/1350 2026-08-18 21:44 by Noways
[教师之家] 为什么余额宝的年化利率越来越低?主要原因有哪些? +3 瞬息宇宙 2026-08-15 3/150 2026-08-18 21:40 by xiaohuixin
[基金申请] 明天放榜? +5 Shxjjxjkx 2026-08-18 5/250 2026-08-18 18:14 by -大大大大大-
[基金申请] 时间戳变了,能看出什么问题? +8 基诺咪客 2026-08-17 11/550 2026-08-18 17:13 by 悠然熊猫
[论文投稿] 投稿咨询 +4 wwm09 2026-08-17 6/300 2026-08-18 15:36 by wwm09
[基金申请] 时间戳又变了8-15 +14 archvillain 2026-08-15 26/1300 2026-08-18 13:37 by phantomgost
[基金申请] 哪位老哥知道今年的国自然具体哪一天放榜? +13 Ldrop2023 2026-08-13 16/800 2026-08-18 12:25 by 淀粉搬运工
[基金申请] 快农历七夕节了,轻松一下,男人悄悄话,女施主请不要进来。 +4 Tide man 2026-08-14 5/250 2026-08-18 11:35 by Tide man
[基金申请] 93BebMhtakh前后11位开头都是大写 +4 且听虎啸 2026-08-17 5/250 2026-08-18 00:49 by 蔡棒棒菂
[基金申请] 感觉是下周放榜了 +6 angus9576 2026-08-17 11/550 2026-08-17 23:57 by angus9576
[基金申请] filecode=后面第一个是大写字母 +8 wangze12014 2026-08-14 10/500 2026-08-17 17:05 by xter9665
[基金申请] 今天系统多次维护,明天很可能放榜! +8 zju2000 2026-08-16 9/450 2026-08-17 12:20 by lmz0216
[精细化工] 招聘 金属平磨液,抛光液研发工程师 +3 小天0311 2026-08-14 3/150 2026-08-16 07:31 by H9PLUS
[基金申请] 各位道友,我要去昆明玩几天,回来见。 +7 Tide man 2026-08-14 8/400 2026-08-15 01:11 by arzu_hma
[基金申请] 重要来源:本周末出结果 +10 瞬息宇宙 2026-08-12 10/500 2026-08-13 15:46 by likettle
[基金申请] 不应该看fileCode +7 且听虎啸 2026-08-12 9/450 2026-08-13 14:27 by flydreamws
[基金申请] Filecode 又变了,巨变 +3 WH3796 2026-08-12 4/200 2026-08-13 14:13 by 小木虫6752397
信息提示
请填处理意见