科学网

 找回密码
  注册

tag 标签: 迭代法

相关帖子

版块 作者 回复/查看 最后发表

没有相关内容

相关日志

Gauss-Seidel 迭代法解线性方程组 A x = b 及其MATLAB编程实现
fzmath 2016-10-23 19:38
Guass-Seidel(高斯--赛德尔) 迭代法(简称 $G-S$ 迭代)是对 Jacobi 迭代的一种改进. 考查 在方程组(6.42) 式, 计算次序是$x_1{(k+1)}\rightarrow x_2^{(k+1)}\rightarrow x_3^{(k+1)}$, 但在计算 $ x_2^{(k+1)}$ 时并没有利用已经计算出的$x_1^{(k+1)}$,而仍利用$x_1^{(k)}$; 同样, 在计算$x_3^{(k+1)}$ 时, 也没有利用最新计算出的 $x_1^{(k+1)}$,$x_2^{(k+1)}$. 为此, 将$(6.42)$ 改进为 这就是G-S迭代. 考虑一般情形, G-S迭代也就是将 J 迭代公式 $(6.40)$ 改进为 这是G-S迭代的分量形式. G-S迭代的矩阵形式即是将 $(6.39)$ 式改为 于是得 即 其中, 迭代矩阵 . Gauss-Seidel 迭代法的MATLAB程序: function = LinGauSeid(A,b,eps,it_max) % 求线性方程组的Gauss-Seidel迭代法,调用格式为 % = LinGauSeid(A,b,ep,it_max) % 其中 % A 为线性方程组的系数矩阵,b 为常数项,ep 为精度要求,默认为1e-5, % it_max 为最大迭代次数,默认为100 % x 为线性方程组的解,k迭代次数 if nargin 4 it_max = 100; end; if nargin 3 eps = 1e-5; end; if min(diag(A))1e-10 error('% 对角元素为0,计算失败!'); end n = length(b); xk = zeros(n,1); k=1; B = -tril(A)\triu(A,1); f = tril(A)\b; while kit_max xk1 = B*xk + f; if max(abs(xk1-xk))eps break; end xk = xk1; k = k+1; end 在命令窗口输入 A = ; b = ; = LinGauSeid(A,b) x = 3.0000 4.0000 -5.0000 k = 25 说G-S迭代是J迭代的改进,其含义是: 当二者均收敛时, G-S迭代比J迭代收敛速度要快. 然而, 也存在J迭代收敛, 而G-S迭代发散的线性方程组.
20697 次阅读|0 个评论
离开计算器可以走多远:x + ln x = 0
热度 1 zjzheng9805 2016-1-19 16:54
离开计算器可以走多远: x + ln x =0 2015 年秋季研究生课程《高等应用数学》期末测验的一道附加题: 求方程 x + ln x = 0 的实根,要求给出分析和计算过程,结果的有效数字越多越好,且不能使用计算器。 在上一年度的期末测验中也有一道类似的估算题,难度稍微要大一些。168 名研究生参加测验,几尽无视附加题,也或许是瞄一眼整个人都不好了,直接放弃。仅有十余名学生动笔写了点东西,而仅有一名学生稍稍走得远一些,但也还是差临门一脚。此次稍稍降低了难度,在确保思想性和技巧性不差的情况下适当地降低了所需手工计算的难度。142 名研究生和2 名本科生参加测验,虽然学生数少了,但约半数学生动了笔,不过都走得不是太远,相比较离目标最近的还是去年的那名学生。考虑到考试时间的限制,而且没有任何提示,能走出一两步也已经很不容易了,但即使时间不限,不提示,不讨论,不会的也还是很可能不会。 本年度首次有本科生选修该研究生课程,自然对这两名本科生的表现较为期待,就谈谈他和她如何处理这道附加题。他先用作图法找到根所在的区间 (0,1) ,然后采用牛顿迭代格式 x k +1 = x k - ( x k + ln x k )/(1 + 1/ x k ), k = 0, 1, 2, ... 给出 x 0 = 1, x 1 = 0.5, x 2 = ? ,然而由于不能使用计算器,他没能走得更远。她给的第一种思路也是牛顿迭代法,并给出了 x 0 = 0.6 和 x 1 = 0.55 ,但未详细说明 0.55 是如何手算出来的。她显然意识到在这个思路下没有计算器是走不了更远的,而且牛顿迭代法是本科掌握的方法,在大学数学复习课上有提到,但不是该课程考察的知识。进而她在第二种思路中试图构造一个小参数 0.1 以采用级数法,然而经过一些尝试发现未能找到合适的构造就放弃了。事实上,在后一种思路中,她已经跨出了大大的一步,若稍加提示,必将天堑变通途。 作为附加题,题目可以较为开放,任何方法都可以尝试,目的也是希望学生能够灵活运用课程中教授的方法。如果不限制使用计算器,牛顿迭代法无疑是可选的最佳方案,迭代效率高。通常的数学软件都采用了牛顿迭代法,如Mathematica 软件的FindRoot 函数或Matlab 软件的fzero 函数。这个时代的人啊,敲几下键盘立马就可以获得精度极高的结果,例如用Mathematica 软件输入FindRoot , { x ,1}] ,可得{ x -0.567143} 。然而,重要的并不是这个附加题的结果。 该课程以林家翘先生和西格尔的《自然科学中确定性问题的应用数学》第一、二卷为主要教材,旨在通过案例讲授“应用数学过程”的思想,贯穿始终介绍了各种各样的摄动方法:级数法、参数微商法、迭代法、庞加莱方法、分部积分法、拉普拉斯方法、傅里叶分析、小波分析、匹配法、WKB 方法、多重尺度法、均匀化方法。解答该附加题也自然可以只限定在这些方法中尝试,容易猜测级数法和迭代法等正则摄动法最为可能可行。 摄动法的思想是将精确解作微小的扰动变成近似解,使得它满足一些可解的方程用以代替原来很难或不能精确求解的方程。因此,对解有初步的了解是很重要的。若把 x + ln x 看作函数,可以采用零值定理找出解所在的大概区间,如 0 x 1 ,当然也可以采用作图法。可是,居然还有更直接的方法,考虑到 x 和 ln x 必须一正一负,立即可得 0 x 1 。这么直接而清晰的一个方法是一名研究生写下的,虽然他写下这点后就没有再走多远。由解所在的区间说明可以将 x 视为小量,那么借助泰勒展开令人讨厌的 ln x 似乎就可以不那么讨厌了。倘若将 ln x 在 x = 0 附近作泰勒展开,即麦克劳林展开,那是行不通的。换个思路,1 - x 也是小量,由泰勒展开(见附A ),可得: ln x = ln ~ - (1 - x ) - (1 - x ) 2 /2 - (1 - x ) 3 /3 - (1 - x ) 4 /4 - (1 - x ) 5 /5 - ..., 至此,可以得到一个能够手工计算的迭代式子,如 x k +1 = 1/2+1/2* , k = 0, 1, 2, ..., 但算起来很辛苦。例如,若保留方括号中两项进行迭代,有: x 0 = 0.5, x 1 = 0.5833, x 2 = 0.5555, x 3 = 0.5640, x 4 = 0.5613, x 5 = 0.5622, ... 辛辛苦苦算来算去也仅能精确到小数点后两位。为提高精度必须保留更高阶的部分,计算难度越发大了去。如果把方程改写为 x *exp( x ) = 1 或 exp( - x ) = x ,将左端作麦克劳林展开再构造迭代式,这样的思路也可行,但也还是会算得很辛苦。至此说明可以通过构造特殊的迭代格式进行计算,与牛顿迭代法相比,效率低一些,但至少是能够手工计算的。实际上,“有效数字越多越好”这个要求的信息量很大,暗示不能蛮干。嗯,还是可以再走得更远一些的。 迭代法通过不断迭代来提高结果的精度,但每次迭代所需的计算量越来越大,而且还有一个严重的缺点,每次迭代都要重复计算出已经精确到的部分。这个缺点是有办法克服的,比方改用级数法,这在课程教学中已经说明过。然而,整个教科书几乎都是围绕着包含小参数的问题在讨论,该附加题所给的方程可不包含小量,那似乎就没辙了。嗯,得先构造一个小量。迭代法可以不用明确小参数,但是可以帮助构造小参数。如果先将方程改写为 x = exp( - x ) ,再运用麦克劳林展开,可得: x = 1 - x + x 2 /2! - x 3 /3! + x 4 /4! - x 5 /5! + ... 然后基于迭代法思想将方程进一步改写为: 2x -1 = x 2 /2! - x 3 /3! + x 4 /4! - x 5 /5! + ... 等式左边的每一项至少是 O (x) ,而右边的每一项都是 o ( x ) 。若略去右边的部分可得 x = 1/2 ,可将其取为小量,记为 a = 1/2 。方程可以进一步改写为: 2 x -2 a = x 2 /2! - x 3 /3! + x 4 /4! - x 5 /5! + ... 进而,可以采用级数法求解得到: x = a + 1/4* a 2 + 1/24* a 3 - 1/192* a 4 - 13/1920* a 5 +... 考虑到 a = 1/2 ,上式的每一项都可以较轻松地通过手工计算完成,各项分别为0.5, 0.0625, 0.0052083, -0.0003255, -0.0002115 等等。因此,若保留前五项可得结果为 0.567171 ,精确到了小数点后四位;项数取得越多精度越高,但收敛速度越来越慢。实际上,也可取小量 c = 1/10 = 0.1 ,正如那名本科生所期待的。进而将方程改写为: 2 x -10 c = x 2 /2! - x 3 /3! + x 4 /4! - x 5 /5! + ... 其级数解中每一项的数量级将非常清晰,见附B 。然而,对于这个题嘛,其实际的计算量没有太多差别,因为收敛速度是一样的。 还能走多远呢?路漫漫其修远兮 ~~~ 郑志军 2016 年1 月7 日 附 A :麦克劳林展开 ln(1 + u ) = u - u 2 /2 + u 3 /3 - u 4 /4 + u 5 /5 - u 6 /6 + ... = Sum exp( u ) = 1 + u + u 2 /2! + u 3 /3! + u 4 /4! + u 5 /5! + u 6 /6! + ... = Sum 附 B :解答小结 1. 方程 x + ln x = 0 的根满足 0 x 1 ,说明 x 可以作为小量 2. 由麦克劳林展开有: x = exp( - x ) = 1 - x + x 2 /2! - x 3 /3! + x 4 /4! - x 5 /5! + ... 3. 由迭代法思想有:2 x -1 = x 2 /2! - x 3 /3! + x 4 /4! - x 5 /5! + ... ,右边的每一项都是 o ( x ) 4. 取小量 c = 1/10 ,将方程进一步改写为:2 x -10 c = x 2 /2! - x 3 /3! + x 4 /4! - x 5 /5! + ... 5. 由级数法可得: x = 5 c + 25/4* c 2 + 125/24* c 3 - 625/192* c 4 - 8125/384* c 5 - 146875/4608* c 6 + 1140625/64512* c 7 + 191171875/1032192* c 8 + ... 6. 手工计算可得: x = 5 c + 6.25 c 2 + 5.2083 c 3 - 3.255 c 4 - 21.15 c 5 - 31.8 c 6 + 17. c 7 + 185. c 8 +... = 0.567143...
个人分类: 应用数学|4839 次阅读|1 个评论
COMSOL求解器 Direct and Iterative Solvers
thomasnust 2013-11-26 12:14
COMSOL求解器详细介绍 推荐- solvers_COMSOL_42.pdf http://www.michelsencentre.com/doc//PDF%20dokumenter/COMSOL/solvers_COMSOL_42.pdf **************************************************** 具体而言- Conditional Number 条件数 Condition number.pdf http://en.wikipedia.org/wiki/Condition_number 直接法 LU分解 LU decomposition.pdf http://en.wikipedia.org/wiki/LU_decomposition MUMPS:a MUltifrontal Massively Parallel sparse direct Solver MUMPS _ a parallel sparse direct solver.pdf http://graal.ens-lyon.fr/MUMPS/ PARDISO:Parallel Direct Sparse Solver Interface PARDISO_ - Parallel Direct Sparse Solver Interface.pdf http://software.intel.com/sites/products/documentation/hpc/mkl/mklman/GUID-7E829836-0FEF-46B2-8943-86A022193462.htm SPOOLES:SParse Object Oriented Linear Equations Solver SPOOLES.pdf http://www.netlib.org/linalg/spooles/spooles.2.2.html 迭代法 共轭梯度 Conjugate gradient Conjugate gradient method.pdf http://en.wikipedia.org/wiki/Conjugate_gradient_method GMRES:Generalized minimal residual method GMRES - Generalized minimal residual method.pdf http://en.wikipedia.org/wiki/Generalized_minimal_residual_method BiCGSTAB:Biconjugate gradient stabilized method BiCGSTAB- Biconjugate gradient stabilized method.pdf http://en.wikipedia.org/wiki/Biconjugate_gradient_stabilized_method
个人分类: 数值方法|12998 次阅读|0 个评论

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-6-3 18:03

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部