| 查看: 3304 | 回复: 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：各位大侠对于这个不动点解法，有没有优化的方法，或者有没有其它可行的解法？谢谢先！ |
回复此楼» 猜你喜欢
基金申请
已经有44人回复
-
CSC与新西兰维多利亚大学PhD奖学金项目
已经有0人回复
-
物理学I论文润色/翻译怎么收费?
已经有128人回复
-
新西兰Robinson研究所 招聘CSC公派访问人员
已经有0人回复
-
帮我的英语口语老师找学生
已经有0人回复
-
什么时候开奖？
已经有13人回复
-
散金币祈福
已经有94人回复
-
青基已中
已经有1人回复
-
散金币祈福
已经有107人回复
-
求助，如何提取ELK的rt-TDDFT在某一时刻的自旋密度分布和电子密度分布
已经有1人回复
【答案】应助回帖
|
这是有约束非线性方程组优化问题，自变量的范围形成一个线性约束。利用内点罚函数算法（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
onesupeng
金虫 (职业作家)
- 计算强帖: 13
- 应助: 256 (大学生)
- 贵宾: 1.36
- 金币: 2116.7
- 散金: 9271
- 红花: 92
- 帖子: 4585
- 在线: 1304.4小时
- 虫号: 394701
- 注册: 2007-06-07
- 专业: 流体力学
4楼2013-12-23 08:47:21
投票:
30