| 查看: 3416 | 回复: 26 | |||
| 当前只显示满足指定条件的回帖,点击这里查看本话题的所有回帖 | |||
[求助]
求助:非线性方程组的求解(郁闷中) 已有10人参与
|
|||
|
各位大侠好,本人工科出身,数学功底实在一般。最近自编一个计算程序,涉及到一个非线性方程组的求解,无奈解法不理想,很多时候不收敛,求大侠指导一二。 方程形式:方程组中的每个方程的形式都是这样的,x + A = f1(x) + f2(x) ,其中f(x)的形式为 f(x) = (x+ B) / ln (x + C) 其中A,B,C为常数。 解法:采用不动点迭代法,即假设一组初值,带入方程的右边,从而得到一组新的值。如误差值大于允许误差,采用加权因子的方式获得新的迭代值,加权因子从0.1到0.625已尝试过多个。 问题:有时迭代过程中变量计算值超过边界条件。例如x的允许范围为8<= x <=12 , 迭代过程中x会超过12或小于8,会导致计算出错。因此我限定如果x超过12,则等于12;小于8,则等于8。但是没有效果,最后x值总是一直超过12就是小于8。费解。 Help:各位大侠对于这个不动点解法,有没有优化的方法,或者有没有其它可行的解法?谢谢先! |
» 猜你喜欢
澳门大学IAPME李海峰教授常年招收博士研究生
已经有0人回复
深圳大学应用技术学院招聘凝聚态物理博士后
已经有124人回复
物理学I论文润色/翻译怎么收费?
已经有149人回复
求助中医院临床中医中标的青年标书学习,谢谢
已经有5人回复
浙大物院表面物理课题组长期招收普博及直博学生
已经有259人回复
【征稿】Springer Nature热电专题 Discover Materials / Discover Chemistry
已经有47人回复
25届硕士申请,材料计算方向
已经有2人回复
课题组招1-2名硕士
已经有1人回复
机时越来越贵之后,我把一部分 DFT 流程搬到了 AI 计算平台上
已经有0人回复
斯坦福大学与Elsevier世界前2%学者新排名查询
已经有0人回复
【答案】应助回帖
|
这是有约束非线性方程组优化问题,自变量的范围形成一个线性约束。利用内点罚函数算法(MATLAB优化工具箱自带)。 以下是如何求数值解。注意,为了简化,让优化工具箱自行计算方程组的数值一阶微分(雅可比矩阵)和二阶微分(Hessian矩阵)。如果结果不理想,可以先调节 options 中的参数。如果还不理想,可以利用符号运算工具箱求出雅可比矩阵和 Hessian 矩阵。 1. 建立一个函数文件 func1.m function y=obj_func(x) % x 为列向量 %常量 A=0.8; B1=0.1; C1=0.4; B2=0.2; C2=0.5; %临时变量 f1 = (x+ B1)./ln(x + C1); f2 = (x+B2)./ln(x+C2); temp_y = x + A - f1 - f2; %目标方程为每项方程的平方和 y = 0.5*sum(temp_y.^2); 2. 建立一个主程序 solve_6_unknowns.m %清屏,清工作区 clc clear all %设置优化算法参数 maxiter = 20000; maxfuneval = length(x0)*maxiter; options = optimset('Display‘, ’off',... ‘Algorithm', 'interior-point',... 'GradObj', 'off',... 'GradConstr', 'off',... 'Hessian', 'off',... 'TolX', 1e-6,... 'TolFun', 1e-6,... 'MaxIter', maxiter,... 'MaxFunEvals', maxfuneval); %设置初始值 x0 = [0.5, 0.5, 0.5, 0.5, 0.5, 0.5]'; %设置边界 lb = [8, 8, 8, 8, 8, 8]'; ub = [12, 12, 12, 12, 12, 12]'; %调用优化函数 [x, fval, exitflag, output] = fmincon(@obj_func, x0, ... [], [], ... [], [], ... lb, ub, ... [], options); %显示结果 disp(['Iterations: ',num2str(output.iterations)]) disp(['Func Evals: ', num2str(output.funcCount)]) disp(['Algorithm: ',output.algorithm]) disp(['exit flag = ',num2str(exitflag)]) disp(['error = ( ',num2str(fval','%-15.6e'),' )']) disp(['x1 = ( ',num2str(x(1),'%-15.6f'),' )']) disp(['x2 = ( ',num2str(x(2),'%-15.6f'),' )']) disp(['x3 = ( ',num2str(x(3),'%-15.6f'),' )']) disp(['x4 = ( ',num2str(x(4),'%-15.6f'),' )']) disp(['x5 = ( ',num2str(x(5),'%-15.6f'),' )']) disp(['x6 = ( ',num2str(x(6),'%-15.6f'),' )']) |
27楼2014-01-08 23:50:34
jerkwin
专家顾问 (正式写手)
-

专家经验: +14 - 计算强帖: 1
- 应助: 454 (硕士)
- 金币: 20659.1
- 散金: 148
- 红花: 81
- 帖子: 813
- 在线: 2648.7小时
- 虫号: 1023452
- 注册: 2010-05-19
- 专业: 理论和计算化学
- 管辖: 分子模拟
2楼2013-12-22 22:11:22
dingd
铁杆木虫 (职业作家)
- 计算强帖: 4
- 应助: 1641 (讲师)
- 金币: 15037.3
- 散金: 101
- 红花: 234
- 帖子: 3410
- 在线: 1223.7小时
- 虫号: 291104
- 注册: 2006-10-28
3楼2013-12-22 22:16:24

4楼2013-12-23 08:47:21










回复此楼
200