科学计算方法解读课件

科学计算的背景非线性方程求根算法线性方程组求解直接法线性方程组求解迭代法《科学计算方法》科学计算的背景《科学计算方法》科学计算方法与计算机有机结合构造出强有力的工作平台数值分析——研究用计算机求解1969年,Apollo

登月计划实现1981年,Columbia号航天飞机发射成功数学问题的方法(算法)和理论方程组求解、方程求根、数据插值、数据拟合、数值积分、微分方程求解vonNeumann1994年,GPS完全投入使用科学计算方法与计算机有机结合数值分析——研究用计算机求解19例1:圆内接正多边形边长计算Pi方法评价算法的主要指标:速度和精度简单迭代算法:

n

Lerror1923.14145241.4e-0043843.14155763.5e-0053.14159264.6e-010例1:圆内接正多边形边长计算Pi方法评价算法的主要指标:速例2.通信卫星覆盖地球面积数学模型实际问题获取数据数值方法、程序数据结果将地球考虑成一个球体,设R为地球半径,h为卫星高度,D为覆盖面在切痕平面上的投影(积分区域)例2.通信卫星覆盖地球面积数学模型实际问题获取数据数值方法

假设某一数据的准确值为

x*,其近似值为x，则称

而称为

x

的相对误差误差的有关概念

e(x)=x

-x*

为

x的绝对误差假设某一数据的准确值为x*,其近似值而称为x的相对如果存在一个适当小的正数ε

，使得

则称ε为绝对误差限。

称εr为相对误差限。

如果存在一个适当小的正数εr

，使得

如果存在一个适当小的正数ε，使得则称ε为绝对误差限。十进制浮点数表示一台微机价格:￥3999.00,

浮点数表示:0.3999×104地球半径:6378137m,(6.378137e+006)

浮点数表示:0.6378137×107光速:2.99792458e+008

浮点数表示:0.299792458×109尾数部阶码部十进制浮点数表示一台微机价格:￥3999.00,尾数部阶码部有效数字概念:取的有限位数如下(≈3.1415926)取

x1=3，误差限不超过0.5;取

x2=3.14,误差限不超过0.005

;若近似值x

的绝对误差限是某一位上的半个单位，该位到x

的第一位非零数字一共有n

位，则称近似值x

有n

位有效数字.

取

x3=3.1416,误差限不超过0.00005

;有效数字概念:取x1=3，误差限不超过0.5;取x2

r

d

例4.水中浮球问题

有一半径r=10cm的球体,密度

=0.638.球体浸入水中后,浸入水中的深度d是多少？

根据阿基米德定律,物体排开水的质量就是水对物体的浮力。整理得:d3–3rd2+4r3

=0例4.水中浮球问题根据阿基米德定律,物体排开水的质量由

=0.638,r=10.代入,得d3–30d2+2552=0令

f(x)=x3–30x2+2552,函数图形如下所示求解方程

f(x)=0,即是求函数

f(x)的零点.f(x)的零点所在区间为:[0,20]roots([1-3002552])ans=26.3146

11.8615-8.1761由=0.638,r=10.代入,得d3–30第一步:对根进行隔离,找出隔根区间,或在隔根区间内确定一个解的近似值x0;设f(x)=0的根为

x*,通过迭代计算,产生序列:

x0

x1

x2

···

xn·········用数值方法求非线性方程的根,分两步进行:第二步:逐步逼近,利用近似解x0(或隔根区间)

通过迭代算法得到更精确的近似解.只须第一步:对根进行隔离,找出隔根区间,或在隔根区间内确定一个解已知方程

f(x)=0有一隔根区间[a,b],且f(x)满足f(a)·f(b)<0,则先将[a,b]等分为两个小区间,判断根属于哪个小区间,舍去无根区间保留有根区间[a1,b1];二分法迭代把区间[a1,b1]一分为二,进一步判断根属于哪个更小的区间[a2,b2],如此不断二分以缩小区间长度.已知方程f(x)=0有一隔根区间[a,b],且f(x)满[a,b]x0=0.5(a+b)[a1,b1]=[a,x0][a1,b1]=[x0,b]x1=0.5(a1+b1)f(a1)f(b1)<0已知f(x)=0在[a,b]内有一根,且f(a)f(b)<0(1)计算:yaf(a),x00.5(a+b),y0f(x0)

判断,若y0=0,则x0是根,否则转下一步;(2)判断,若y0·ya<0,则a1a,b1

x0

否则

a1x0,b1b,yay0[a,b]x0=0.5(a+b)[a1,b1]=[a,x0二分法迭代将得到一系列隔根区间

定理2.2

设x*是

f(x)=0在[a,b]内的唯一根,且

f(a)·f(b)<0,则二分计算过程中,各区间的中点数列性质:1.f(an)·f(bn)<0;2.bn–

an=(b–a)/2n满足:|xn–x*|≤(b–a)/2n+1思考:

|xn+1–xn

|=？二分法迭代将得到一系列隔根区间定理2.2设x*是f(构造有效的迭代格式选取合适的迭代初值对迭代格式进行收敛性分析例5.圆周率计算初值:

x0=1(n=1,2,3,······)迭代格式:将一个计算过程反复进行称为迭代,迭代法是一类常见常用的计算技术构造有效的迭代格式例5.圆周率计算初值:x0=1(n=初值:x1=1.5迭代格式:xn+1=0.5(xn+2/xn)(n=1,2,·····)例6.平方根算法求xn

Error1.4166666666666672.45e-0031.4142156862745102.12e-0061.4142135623746901.59e-0121.4142135623730952.22e-0161.4142135623730952.22e-016初值:x1=1.5例6.平方根算法求设

x*是方程

f(x)=0的根,x0是x*的近似值.在

x0附近,对函数做局部线性化x1比x0更接近于x*x0x1x*f(x)=0设x*是方程f(x)=0的根,x0是x*的近似(n=0,1,2,·····)牛顿迭代格式给定初值

x0,迭代产生数列x0,x1,x2,·········,

xn,

·······应用——求正数平方根算法设C>0,x2–C=0令

f(x)=x2–C,则(n=0,1,2,·····)牛顿迭代格式给定初值由此可知,平方根迭代具有

2阶收敛速度

由此可知,平方根迭代具有2阶收敛速度f’<0,f”>0f’>0,f”>0f’>0,f”<0f’<0,f”<0牛顿迭代法收敛的四种情况f’<0,f”>0f’>0,f”>0f’>0,例7.牛顿迭代法的收敛域问题:

用牛顿迭代法求解复数方程

z3–1=0，该方程在复平面上三个根分别是z1=1选择中心位于坐标原点，边长为2的正方形内的任意点作初始值，进行迭代，把收敛到三个根的初值分为三类，并分别标上不同颜色（例如红、黄、蓝）。对充分多的初始点进行实验，绘出牛顿迭代法对该方程的收敛域彩色图。

例7.牛顿迭代法的收敛域问题:z1=1选择中心位于坐标收敛到z1的牛顿迭代初值点集合收敛到z2的牛顿迭代初值点集合收敛到z3的牛顿迭代初值点集合收敛到z1的牛顿迭代初值点集合收敛到z2的牛顿迭代初在复平面内,有一些例外点是牛顿迭代不收敛的初值点.这些例外点构成了茹利亚集(为纪念法国女数学家Julia).在复平面内,有一些例外点是牛顿迭代不收敛的初值点.这些例外线性方程组的矩阵形式a11x1+a12x2+····+a1nxn=b1a21x1+a22x2+····+a2nxn=b2·········································an1x1+an2x2+····+annxn=bnAX=b(i=1,2,···,n)线性方程组求解:1.直接方法;2.基本迭代法;3.子空间方法X

?

b线性方程组的矩阵形式a11x1+a12x2+····+a例8x4=2,x3=4,x2=–1,x1=3例8x4=2,x3=4,x2=–1,x1=3解上三角方程组计算：xn

=bn

/ann

(a11…ann≠0)xk=[bk－(ak,k+1xk+1+…+ak

n)]/akk

(k=n－1,···,1)除法:n次;乘法:n(n-1)/2次,乘、除法运算共n(n+1)/2次,简记为

O(n2)解上三角方程组计算：xn=bn/ann(a11…an消元过程(化一般方程组为上三角方程组)消元过程(化一般方程组为上三角方程组)增广矩阵

计算:[m21

m31

m41]T=[a21

a31

a41]T/a11

用–m21乘矩阵第一行加到矩阵第二行;用–m31乘矩阵第一行加到矩阵第三行;用–m41乘矩阵第一行加到矩阵第四行;增广矩阵计算:[m21m31m41]T=[a实现第一轮消元实现第二轮消元、第三轮消元·········消元过程乘除法次数:O(n3)实现第一轮消元实现第二轮消元、第三轮消元·········消例9.特点:系数矩阵主对角元均不为零计算格式

X(1)=BX(0)+f取

X(0)=例9.特点:系数矩阵主对角元均不为零计算格式X(1)=

X*1.00001.00001.0000

X(0)000

X(1)

0.77780.80000.8667

X(2)0.96300.96440.9778

X(3)0.99290.99350.9952计算格式:X(k+1)=BX(k)+fX(0),X(1),···X(n)·····

X(4)········0.99870.99880.9991X(n+1)≈

BX(n)+fX*X(0)X(1)X(2)X雅可比迭代法(i=1,2,…n;k=1,2,……)取初始向量X(0)=[x1(0)x2(0)···xn(0)]T,迭代计算(i=1,2,…,n)雅可比迭代法(i=1,2,…n;k=1,2,……)取初雅可比迭代法的矩阵表示将方程组AX=b

的系数矩阵

A

分解

A=D–U–LAX=b=>DX(k+1)=(U+L)X(k)+bX(k+1)=D-1(U+L)X(k)+D-1b记BJ=D-1(U+L)X(k+1)=BJX(k)+fJ雅可比迭代法的矩阵表示将方程组AX=b的系数矩阵A雅可比迭代矩阵雅可比迭代矩阵例10误差限5e-0045e-0055e-0065e-007赛德尔迭代次数5567雅可比迭代次数67910高斯-赛德尔迭代法例10误差限5e-0045e-0055e-0065e-007高斯-赛德尔迭代法矩阵表示AX=bA=M–NMX(k+1)=NX(k)+b高斯-赛德尔迭代法矩阵表示AX=bA=M–(i=1,2,…,n)高斯-赛德尔迭代法计算格式(i=1,2,…n;k=1,2,……)取初始向量x(0)=[x1(0)x2(0)···xn(0)]T,

做迭代计算(i=1,2,…,n)高斯-赛德尔迭代法计算格式(i=AX=b(M–N)X=b记

(k)=X(k)–X*(k=0,1,2,3,······)则有

(k+1)=B(k)(k=0,1,2,3,······)计算格式:X(k+1)=BX(k)+f(B=M-1N)

X(k+1)–X*=B(X(k)–X*)设方程组的精确解为X*,则有X*=BX*+fAX=b(M–N)X=b记平面点列:Xk∈Rn:X1,X2,···,Xk,···················平面点列:Xk∈Rn:X1,X2,···,证:由(k)=B(k-1),得

||(k)||≤||B||||(k-1)||

(k=1,2,3,······)所以命题

若||B||<1,则迭代法

X(k+1)=BX(k)+f

收敛||(k)||≤||B||k||(0)||

||B||<1证:由(k)=B(k-1),得所以命题若|定理4.3

若Ax=b的系数矩阵A是严格对角占优矩阵,则Jacobi迭代和Seidel迭代均收敛定理4.2:设X*为方程组

AX=b的解若||B||<1,则对迭代格式

X(k+1)=BX(k)+f

有(1)(2)定理4.3若Ax=b的系数矩阵A是严格对角占优矩阵,则Ja例11.平面温度场问题:令

h=1/(n+1),xj=jh,yj=jh(i,j=0,1,···,n+1)记

ui,j=u(xi,yj),(i,j=0,1,···,n+1)迭代格式(i,j=1,···,n)u0,j=0,ui,0=0,ui,n+1

=0例11.平面温度场问题:令h=1/(n+1),结点数n2102202

402迭代次数

1826062077CPU时间(s)0.974.32858.531误差

0.00236.4274e-41.6814e-4高斯-赛德尔迭代法实验(误差限10-8):结点数n2102(i=1,2,···,n;k=1,2,3,··········)迭代格式超松驰迭代方法SOR(successiveoverrelaxation)Prof.DavidM.Young1954美国数学科学学报(i=1,2,···,n;k=1,2,3,···Seidel迭代格式SOR迭代格式最佳松驰因子平面温度场的计算问题Seidel迭代格式SOR迭代格式最佳松驰因子平面温度场的计结点数n2102202

402迭代次数

1826062077CPU时间(s)0.974.32858.531误差

0.00236.4274e-41.6814e-4高斯-赛德尔迭代(误差限10-8):SOR迭代实验(误差限10-8):结点数n2102202

402迭代次数

4074137CPU时间(s)0.110.65604.9530误差

0.00236.4306e-41.6944e-4结点数n2102人有了知识，就会具备各种分析能力，明辨是非的能力。所以我们要勤恳读书，广泛阅读，古人说“书中自有黄金屋。”通过阅读科技书籍，我们能丰富知识，培养逻辑思维能力；通过阅读文学作品，我们能提高文学鉴赏水平，培养文学情趣；通过阅读报刊，我们能增长见识，扩大自己的知识面。有许多书籍还能培养我们的道德情操，给我们巨大的精神力量，鼓舞我们前进。人有了知识，就会具备各种分析能力，科学计算方法解读课件科学计算的背景非线性方程求根算法线性方程组求解直接法线性方程组求解迭代法《科学计算方法》科学计算的背景《科学计算方法》科学计算方法与计算机有机结合构造出强有力的工作平台数值分析——研究用计算机求解1969年,Apollo

登月计划实现1981年,Columbia号航天飞机发射成功数学问题的方法(算法)和理论方程组求解、方程求根、数据插值、数据拟合、数值积分、微分方程求解vonNeumann1994年,GPS完全投入使用科学计算方法与计算机有机结合数值分析——研究用计算机求解19例1:圆内接正多边形边长计算Pi方法评价算法的主要指标:速度和精度简单迭代算法:

n

Lerror1923.14145241.4e-0043843.14155763.5e-0053.14159264.6e-010例1:圆内接正多边形边长计算Pi方法评价算法的主要指标:速例2.通信卫星覆盖地球面积数学模型实际问题获取数据数值方法、程序数据结果将地球考虑成一个球体,设R为地球半径,h为卫星高度,D为覆盖面在切痕平面上的投影(积分区域)例2.通信卫星覆盖地球面积数学模型实际问题获取数据数值方法

假设某一数据的准确值为

x*,其近似值为x，则称

而称为

x

的相对误差误差的有关概念

e(x)=x

-x*

为

x的绝对误差假设某一数据的准确值为x*,其近似值而称为x的相对如果存在一个适当小的正数ε

，使得

则称ε为绝对误差限。

称εr为相对误差限。

如果存在一个适当小的正数εr

，使得

如果存在一个适当小的正数ε，使得则称ε为绝对误差限。十进制浮点数表示一台微机价格:￥3999.00,

浮点数表示:0.3999×104地球半径:6378137m,(6.378137e+006)

浮点数表示:0.6378137×107光速:2.99792458e+008

浮点数表示:0.299792458×109尾数部阶码部十进制浮点数表示一台微机价格:￥3999.00,尾数部阶码部有效数字概念:取的有限位数如下(≈3.1415926)取

x1=3，误差限不超过0.5;取

x2=3.14,误差限不超过0.005

;若近似值x

的绝对误差限是某一位上的半个单位，该位到x

的第一位非零数字一共有n

位，则称近似值x

有n

位有效数字.

取

x3=3.1416,误差限不超过0.00005

;有效数字概念:取x1=3，误差限不超过0.5;取x2

r

d

例4.水中浮球问题

有一半径r=10cm的球体,密度

=0.638.球体浸入水中后,浸入水中的深度d是多少？

根据阿基米德定律,物体排开水的质量就是水对物体的浮力。整理得:d3–3rd2+4r3

=0例4.水中浮球问题根据阿基米德定律,物体排开水的质量由

=0.638,r=10.代入,得d3–30d2+2552=0令

f(x)=x3–30x2+2552,函数图形如下所示求解方程

f(x)=0,即是求函数

f(x)的零点.f(x)的零点所在区间为:[0,20]roots([1-3002552])ans=26.3146

11.8615-8.1761由=0.638,r=10.代入,得d3–30第一步:对根进行隔离,找出隔根区间,或在隔根区间内确定一个解的近似值x0;设f(x)=0的根为

x*,通过迭代计算,产生序列:

x0

x1

x2

···

xn·········用数值方法求非线性方程的根,分两步进行:第二步:逐步逼近,利用近似解x0(或隔根区间)

通过迭代算法得到更精确的近似解.只须第一步:对根进行隔离,找出隔根区间,或在隔根区间内确定一个解已知方程

f(x)=0有一隔根区间[a,b],且f(x)满足f(a)·f(b)<0,则先将[a,b]等分为两个小区间,判断根属于哪个小区间,舍去无根区间保留有根区间[a1,b1];二分法迭代把区间[a1,b1]一分为二,进一步判断根属于哪个更小的区间[a2,b2],如此不断二分以缩小区间长度.已知方程f(x)=0有一隔根区间[a,b],且f(x)满[a,b]x0=0.5(a+b)[a1,b1]=[a,x0][a1,b1]=[x0,b]x1=0.5(a1+b1)f(a1)f(b1)<0已知f(x)=0在[a,b]内有一根,且f(a)f(b)<0(1)计算:yaf(a),x00.5(a+b),y0f(x0)

判断,若y0=0,则x0是根,否则转下一步;(2)判断,若y0·ya<0,则a1a,b1

x0

否则

a1x0,b1b,yay0[a,b]x0=0.5(a+b)[a1,b1]=[a,x0二分法迭代将得到一系列隔根区间

定理2.2

设x*是

f(x)=0在[a,b]内的唯一根,且

f(a)·f(b)<0,则二分计算过程中,各区间的中点数列性质:1.f(an)·f(bn)<0;2.bn–

an=(b–a)/2n满足:|xn–x*|≤(b–a)/2n+1思考:

|xn+1–xn

|=？二分法迭代将得到一系列隔根区间定理2.2设x*是f(构造有效的迭代格式选取合适的迭代初值对迭代格式进行收敛性分析例5.圆周率计算初值:

x0=1(n=1,2,3,······)迭代格式:将一个计算过程反复进行称为迭代,迭代法是一类常见常用的计算技术构造有效的迭代格式例5.圆周率计算初值:x0=1(n=初值:x1=1.5迭代格式:xn+1=0.5(xn+2/xn)(n=1,2,·····)例6.平方根算法求xn

Error1.4166666666666672.45e-0031.4142156862745102.12e-0061.4142135623746901.59e-0121.4142135623730952.22e-0161.4142135623730952.22e-016初值:x1=1.5例6.平方根算法求设

x*是方程

f(x)=0的根,x0是x*的近似值.在

x0附近,对函数做局部线性化x1比x0更接近于x*x0x1x*f(x)=0设x*是方程f(x)=0的根,x0是x*的近似(n=0,1,2,·····)牛顿迭代格式给定初值

x0,迭代产生数列x0,x1,x2,·········,

xn,

·······应用——求正数平方根算法设C>0,x2–C=0令

f(x)=x2–C,则(n=0,1,2,·····)牛顿迭代格式给定初值由此可知,平方根迭代具有

2阶收敛速度

由此可知,平方根迭代具有2阶收敛速度f’<0,f”>0f’>0,f”>0f’>0,f”<0f’<0,f”<0牛顿迭代法收敛的四种情况f’<0,f”>0f’>0,f”>0f’>0,例7.牛顿迭代法的收敛域问题:

用牛顿迭代法求解复数方程

z3–1=0，该方程在复平面上三个根分别是z1=1选择中心位于坐标原点，边长为2的正方形内的任意点作初始值，进行迭代，把收敛到三个根的初值分为三类，并分别标上不同颜色（例如红、黄、蓝）。对充分多的初始点进行实验，绘出牛顿迭代法对该方程的收敛域彩色图。

例7.牛顿迭代法的收敛域问题:z1=1选择中心位于坐标收敛到z1的牛顿迭代初值点集合收敛到z2的牛顿迭代初值点集合收敛到z3的牛顿迭代初值点集合收敛到z1的牛顿迭代初值点集合收敛到z2的牛顿迭代初在复平面内,有一些例外点是牛顿迭代不收敛的初值点.这些例外点构成了茹利亚集(为纪念法国女数学家Julia).在复平面内,有一些例外点是牛顿迭代不收敛的初值点.这些例外线性方程组的矩阵形式a11x1+a12x2+····+a1nxn=b1a21x1+a22x2+····+a2nxn=b2·········································an1x1+an2x2+····+annxn=bnAX=b(i=1,2,···,n)线性方程组求解:1.直接方法;2.基本迭代法;3.子空间方法X

?

b线性方程组的矩阵形式a11x1+a12x2+····+a例8x4=2,x3=4,x2=–1,x1=3例8x4=2,x3=4,x2=–1,x1=3解上三角方程组计算：xn

=bn

/ann

(a11…ann≠0)xk=[bk－(ak,k+1xk+1+…+ak

n)]/akk

(k=n－1,···,1)除法:n次;乘法:n(n-1)/2次,乘、除法运算共n(n+1)/2次,简记为

O(n2)解上三角方程组计算：xn=bn/ann(a11…an消元过程(化一般方程组为上三角方程组)消元过程(化一般方程组为上三角方程组)增广矩阵

计算:[m21

m31

m41]T=[a21

a31

a41]T/a11

用–m21乘矩阵第一行加到矩阵第二行;用–m31乘矩阵第一行加到矩阵第三行;用–m41乘矩阵第一行加到矩阵第四行;增广矩阵计算:[m21m31m41]T=[a实现第一轮消元实现第二轮消元、第三轮消元·········消元过程乘除法次数:O(n3)实现第一轮消元实现第二轮消元、第三轮消元·········消例9.特点:系数矩阵主对角元均不为零计算格式

X(1)=BX(0)+f取

X(0)=例9.特点:系数矩阵主对角元均不为零计算格式X(1)=

X*1.00001.00001.0000

X(0)000

X(1)

0.77780.80000.8667

X(2)0.96300.96440.9778

X(3)0.99290.99350.9952计算格式:X(k+1)=BX(k)+fX(0),X(1),···X(n)·····

X(4)········0.99870.99880.9991X(n+1)≈

BX(n)+fX*X(0)X(1)X(2)X雅可比迭代法(i=1,2,…n;k=1,2,……)取初始向量X(0)=[x1(0)x2(0)···xn(0)]T,迭代计算(i=1,2,…,n)雅可比迭代法(i=1,2,…n;k=1,2,……)取初雅可比迭代法的矩阵表示将方程组AX=b

的系数矩阵

A

分解

A=D–U–LAX=b=>DX(k+1)=(U+L)X(k)+bX(k+1)=D-1(U+L)X(k)+D-1b记BJ=D-1(U+L)X(k+1)=BJX(k)+fJ雅可比迭代法的矩阵表示将方程组AX=b的系数矩阵A雅可比迭代矩阵雅可比迭代矩阵例10误差限5e-0045e-0055e-0065e-007赛德尔迭代次数5567雅可比迭代次数67910高斯-赛德尔迭代法例10误差限5e-0045e-0055e-0065e-007高斯-赛德尔迭代法矩阵表示AX=bA=M–NMX(k+1)=NX(k)+b高斯-赛德尔迭代法矩阵表示AX=bA=M–(i=1,2,…,n)高斯-赛德尔迭代法计算格式(i=1,2,…n;k=1,2,……)取初始向量x(0)=[x1(0)x2(0)···xn(0)]T,

做迭代计算(i=1,2,…,n)高斯-赛德尔迭代法计算格式(i=AX=b(M–N)X=b记

(k)=X(k)–X*(k=0,1,2,3,······)则有

(k+1)=B(k)(k=0,1,2,3,······)计算格式:X(k+1)=BX(k)+f(B=M-1N)

X(k+1)–X*=B(X(k)–X*)设方程组的精确解为X*,则有X*=BX*+fAX=b(M–N)X=b记平面点列:Xk∈Rn:X1,X2,···,Xk,···················平面点列:Xk∈Rn:X1,X2,···,证:由(k)=B(k-1),得

||(k)||≤||B||||(k-1)||

(k=1,2,3,······)所以命题

若||B||<1,则迭代法

X(k+1)=BX(k)+f

收敛||(k)||≤||B||k||(0)||