第4章数值计算

合集下载

数值分析-第4章数值积分和数值微分

A0+A1=2 A0x0+A1x1=0 A0x02+A1x12=2/3 A0x03+A1x13=0
A0 A1 1 解得: 1 x 0 x1 3
求积公式为
1 1 1 f ( x)dx f ( ) f ( ) 3 3
x f(x)
数值分析
1 4
2 4.5
3 6
4 8
5 8.5
1
一、数值积分的基本概念求积节点数值积分定义如下：是离散点上的函数值的线性组合
I [ f ] f ( x)dx I n [ f ] Ai f ( xi )
b a i 0 n
称为数值积分公式
称为求积系数,与f (x)无关,与积分区间和求积节点有关
b a
Rn ( x) dx
定理：形如 Ak f ( xk ) 的求积公式至少有 n 次代数精度
A 该公式为插值型（即： k a l k ( x)dx ）
数值分析
b
5
例1 试确定参数A0,A1,A2,使求积公式
1 f ( x)dx A0 f (1) A1 f (0) A2 f (1)
证明因为Simpson公式对不高于三次的多项式精确成立。即

b
a
p 2 ( x)dx
ba ab [ p 2 (a) 4 p 2 ( ) p 2 (b)] 6 2
构造三次多项式H3(x),使满足 H3(a)=(a) ,H3(b)=(b),
H 3 (( a b) / 2) f (( a b) / 2), H 3 (( a b) / 2) f (( a b) / 2), 这时插值误差为
1

MATLAB-第4章

v
i 1
n
2 i
。

max { vi } 。
1 ≤i ≤n
设 A 是一个 m ×n 的矩阵，矩阵的 3 种常用范数如下。 1-范数： A 1 max { aij } 。
1 ≤ j ≤n i 1 m
2-范数： A 2 1 ，其中 λ 1 为 A'A 最大特征值。 ∞-范数： A max { aij } 。
【例4.6】先建立5 × 5矩阵A，然后将A的第一行元素乘以1，第二行乘以2，…，第五行乘以5。用一个对角矩阵左乘一个矩阵时，相当于用对角阵的第一个元素乘以该矩阵的第一行，用对角阵的第二个元素乘以该矩阵的第二行……依此类推，因此，只需按要求构造一个对角矩阵D，并用D左乘A即可。命令如下： A=[1:5;2:6;3:7;4:8;5:9] D=diag(1:5); D*A %用D左乘A,对A的每行乘以一个指定常数
（2）构造对角矩阵设V为具有m个元素的向量，diag(V,k)的功能是产生一个 n × n（n = m + k|）对角阵，其第k条对角线的元素即为向量V的元素。例如： diag(1:3,-1) ans = 0 0 0 0 1 0 0 0 0 2 0 0 0 0 3 0 省略k时，相当于k为0，其主对角线元素即为向量V的元素。
2．矩阵的秩与迹（1）矩阵的秩 rank(A) （2）矩阵的迹矩阵的迹即矩阵的对角线元素之和。 trace(A)。
3．向量和矩阵的范数
设向量 V = (v1 ,v2 ,…,vn )，向量的 3 种常用范数如下。 1-范数： V 2-范数： V ? -范数： V
1
vi 。
i 1
n
2

3．矩阵的转置所谓转置，即把源矩阵的第一行变成目标矩阵第一列，第二行变成第二列……依此类推。显然，一个m行n列的矩阵经过转置运算后，变成一个n行m列的矩阵。MATLAB中，转置运算符是单撇号(')。

第4章数值积分与数值微分

1 (a b).得到的求积公式就是中矩形公式。再令 2 f ( x) x 2 , 代入(1.4)式的第三式有
b ab 2 ba 2 1 2 A x (b a)( ) (a b ) x 2 dx (b 3 a 3 ), a 2 4 3 说明中矩形公式对 ( x) x 2不精确成立，故它的代 f 数精确度为 . 1
定理 1 求积公式 f ( x)dx Ak f k 至少具有n次代数精度
a k 0
它是插值型求积公式 .
四、求积公式的余项若求积公式(1.3)的代数精确度为m，则由求积公式余项的表达式(1.7)可以证明余项形如
R[ f ] f ( x)dx Ak f ( xk ) Kf ( m1) ( ), (1.8)
k 0 n
Hale Waihona Puke 第4章数值积分与数值微分
~ 定义 3 若 0, 0,只要 f ( xk ) f k (k 0,, n), 就有 ~ | I n ( f ) I n ( f ) |
《数值分析》
~ Ak [ f ( xk ) f ( xk )] ,
此求积公式的余项。
第4章数值积分与数值微分
1 A1 B0 2 1 1 《 A1 0 x 2 dx 3 1 2 数 1 1 A1 , A0 , B0 于是有 f ( x)dx 2 f (0) 1 f (1) 1 f ' (0) 值解得 3 3 6 3 3 6 分 0 1 1 析当 3时 x 3 dx . 而上式右端为1/3，故公式对 f ( x) x 》 4 0
k 0
n
则称求积公式 (1.3) 是稳定的 .

数值分析第四章外推法计算数值微分MATLAB计算实验报告

数值分析第四章外推法计算数值微分MATLAB计算实验报告数值分析MATLAB计算实验报告姓名班级学号⼀、实验名称⽤MATLAB编程实现数值微分的外推法计算。

⼆、实验⽬的1．掌握数值微分和定义和外推法的计算过程；2．了解数值微分外推法的计算⽅法并且编写出与其算法对应的MATLAB程序代码；3．体会利⽤MATLAB软件进⾏数值计算。

三、实验内容⽤外推法计算f(x)=x2e?x在x=0.5的导数。

四、算法描述1.命名函数。

2.如果输⼊未知数少于四个，默认精度10^-33.描述T表矩阵坐标4.依次赋值计算 T表第⼀列5.根据数值微分计算公式求出T表矩阵的值6.若达到精度则运算结束，若未达到循环计算7.输出T表，得出的值就是导数值五、实验结果六、实验结果分析此实验通过MATLAB实现外推法数值微分计算，得到相应的数据，⽅便对数据进⾏分析。

从结果可以看出，当步长h=0.025时⽤中点微分公式只有3位有效数字，外推⼀次达到5位有效数字，外推两次达到9位有效数字。

七、附录（程序）function g=waituifa(fname,x,h,e)if nargin<4,e=1e-3;end;i=1;j=1;G(1,1)=(feval(fname,x+h)-feval(fname,x-h))/(2*h);G(i+1,1)=(feval(fname,x+h/2)-feval(fname,x-h/2))/h;G(i+1,j+1)=(4^j*G(i+1,j)-G(i,j))/(4^j-1);while abs(G(i+1,i+1)-G(i+1,i))>ei=i+1;G(i+1,1)=(feval(fname,x+h/2^i)-feval(fname,x-h/2^i))/(2*h/2^i); for j=1:iG(i+1,j+1)=((4^j)*G(i+1,j)-G(i,j))/(4^j-1);endendGg=G(i+1,i+1);。

第4章导热问题的数值解法共30页

若取上面式右边的前三项，并将式①和式③相加移项整理即得二阶导数的中心差分：
2t tm 1 ,n2 tm ,ntm 1 ,no( x2)
x2m ,n
x2
截断误差
同样可得：
未明确写出的级数余项中
的Δx的最低阶数为2
2t tm ,n 12tm ,ntm ,n 1o( y2)
y2m ,n
y2
28.05.2020 - 8 -
(3) 实验法: 是传热学的基本研究方法，a 适应性不好； b 费用昂贵
数值解法：有限差分法（finite-difference）、有限元法（finite-element）、边界元法（boundary- element）、分子动力学模拟（MD）
28.05.2020 - 2 -
第4章导热问题的数值解法——§4-1 导热问题数值求解的基本思
2 例题条件
y
h3t f
W
t0
t2 x 2
t2 y 2
0
x 0, t t0
x H,
t x
h2 (t
tf)
h2t f
y 0,
t y
h1 (t
tf)
yW ,
t y
h3 (t
tf)
h1t f
Hx
二维矩形域内稳态无内热源，
常物性的导热问题
28.05.2020 - 4 -
第4章导热问题的数值解法———§4-1 导热问题数值求解的基本思想
第4章导热问题的数值解法———§4-1 导热问题数值求解的基本思想
以二维、稳态、有内热源的导热问题为例此时：
Φ 上 Φ 下 Φ 左＋ Φ 右 Φ v 0 左Ad dxtyd dxt

Ch04：数值计算方法之函数增量的计算方法

6． sin(u)/u的数学解释
对于适当的u，注意到sin(u)/u可看成是sin(x)在x=0处的 u增量比
sin( u ) sin( 0 u ) sin( 0) u u
所以我们采用统一的记号把它记为 sin( u ) DQSINTV (u ) u 利用泰勒展式求数值解的源代码见教材第90页程序4.03。对任何一个连续可微的实函数f(x)，如果我们得到了它在x=x0处的泰勒展式，那么我们不难直接利用泰勒展式求f(x)在x=x0 处的Δx增量比。
1. 问题的提出
当|u|充分小时，对于计算ln (1+u)的数字结果这个数学问题来说，由于数学问题的条件数会变得充分大，所以u的较小的相对误差会导致ln(1+u)的较大的相对误差。利用流行的程序设计语言提供的库函数来计算，不难验证，这个问题依然普遍存在。利用第3章程序3.14来计算，令a=1+u，那么程序中计算x=(a-1)/(a+1)也是非常不利的运算。所以，如果有效地解决了当|u|充分小时ln (1+u)的精确计算问题，那么第3章的遗留问题和我们这一节所要解决的问题都得到解决。
4.4 对数函数增量的计算方法
记LOG(x)为数学中的以e为底的对数函数ln x，利用我们习惯的记号，LOG(x)在x0处的Δx增量和相应的增量比可以记为 DLOG(x0,Δx) =LOG(x0+Δx)-LOG(x0) DQLOG(x0,Δx)=[LOG (x0+Δx)- LOG(x0)]/Δx 虽然利用对数函数的性质可以得到 DLOG(x0,Δx)=LOG (1+Δx/ x0) DQLOG(x0,Δx)=LOG (1+Δx/ x0)/Δx 在许多情况下也能得到非常理想的结果，但还是存在一些问题。

【推荐】数值计算方法：第4章-多项式插值方法.ppt

两点
多项式插值就是直线
，经过这两点的
称给定
为线性插值多项式。称
为关于点
的线性插值基函数，其在节点处满足：
6
4.2.1 线性插值与二次插值假定插值节点为，，，要求二次插值多项式
几何上
是通过三点
可以用基函数的方法求的表源自式，是二次函数，的抛物线.
7
4.2.2 拉格朗日插值多项式
求n+1个次数满足
且次数不超过n 的多项式，其所给出形式的系数为
称
为牛顿（Newton）均差插值多项式.
系数就是均差表4-1中主对角线上的各阶均差，它比拉格朗日插值计算量省，且便于程序设计.
25
4.3.2 Newton均差插值多项式（*）为插值余项，由插值多项式惟一性知，它与
拉格朗日插值多项式的余项应该是等价的. 事实上，利用均差与导数关系式就可以证明这一点. 但（3.7）更有一般性，它在是由离散点（给3.出7）的
式求 x 的近似值。
解 (1) 选取节点x=2,3,4
xf 一二三
kk
(x k)
阶均
阶均
阶均
31
32
4.4 分段低次插值
4.4.1 Runge现象在次数增加时逼近的精度是否也增加?
问题：根据区间上给出的节点做出的插值多项式
事实上，对于有些函数，插值多项式次数很高时会在某些区间内产生较大的误差。例如著名的Runge现象。
分段插值的基本思想是将插值区间划分为若干个小区间, 然后在每个小区间上做满足一定条件的低阶插值.
35
4.4.2 分段低次插值
例如分段线性插值。所谓分段线性插值就是通过插值点用折线段连接起来

数值分析--第4章数值积分与数值微分[1]详解

第4章数值积分与数值微分1 数值积分的基本概念实际问题当中常常需要计算定积分。

在微积分中，我们熟知，牛顿-莱布尼兹公式是计算定积分的一种有效工具,在理论和实际计算上有很大作用。

对定积分()ba I f x dx =⎰，若()f x 在区间[,]ab 上连续，且()f x 的原函数为()F x ,则可计算定积分()()()baf x dx F b F a =-⎰似乎问题已经解决，其实不然。

如1）()f x 是由测量或数值计算给出数据表时,Newton-Leibnitz 公式无法应用。

2)许多形式上很简单的函数，例如222sin 1(),sin ,cos ,,ln x x f x x x e x x-= 等等，它们的原函数不能用初等函数的有限形式表示。

3）即使有些被积函数的原函数能通过初等函数的有限形式表示，但应用牛顿—莱布尼兹公式计算，仍涉及大量的数值计算，还不如应用数值积分的方法来得方便，既节省工作量，又满足精度的要求。

例如下列积分241arc 1)arc 1)1dx tg tg C x ⎡⎤=+++-+⎣⎦+⎰ 对于上述这些情况,都要求建立定积分的近似计算方法—-数值积分法。

1。

1 数值求积分的基本思想根据以上所述,数值求积公式应该避免用原函数表示，而由被积函数的值决定.由积分中值定理：对()[,]f x C a b ∈，存在[,]a b ξ∈，有()()()baf x dx b a f ξ=-⎰表明,定积分所表示的曲边梯形的面积等于底为b a -而高为()f ξ的矩形面积（图4-1）。

问题在于点ξ的具体位置一般是不知道的，因而难以准确算出()f ξ。

我们将()f ξ称为区间[,]a b 上的平均高度。

这样，只要对平均高度()f ξ提供一种算法，相应地便获得一种数值求积分方法.如果我们用两端的算术平均作为平均高度()f ξ的近似值，这样导出的求积公式[()()]2b aT f a f b -=+ (4—1）便是我们所熟悉的梯形公式（图4-2）。

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

例4.3 求积分 0 e dx 。先建立一个函数文件ex.m： function ex=ex(x) ex=exp(-x.^2); %注意应用点运算 return 然后，在MATLAB命令窗口，输入命令： quad(‘ex’,0,1,1e-6) %注意函数名应加字符引号
2
1
(2)被积函数由一个表格定义 MATLAB中，对由表格形式定义的函数关系的求定积分问题用trapz(X,Y)函数。其中向量X、Y定义函数关系Y=f(X)。例4.4 用trapz函数计算积分。在MATLAB命令窗口，输入命令： X=0:0.01:1;Y=exp(-X.^2); trapz(X,Y)
4.4 多项式计算
1. 多项式的建立已知一个多项式的全部根X求多项式系数的函数是 poly(X)，该函数返回以X为全部根的一个多项式P，当X是一个长度为m的向量时，P是一个长度为 m+1的向量。 2. 多项式求根求多项式p(x)的根的函数是roots(P)，这里，P是 p(x)的系数向量，该函数返回方程p(x)=0的全部根 (含重根，复根)。 3. 多项式求值求多项式p(x)在某点或某些点的函数值的函数是 polyval(P,x)。若x为一数值，则求多项式在该点的值；若x为向量或矩阵，则对向量或矩阵中的每个元素求其多项式的值。
2）矩阵的范数及其计算函数
MATLAB中提供了求3种矩阵范数的函数，其函数调用格式与求向量的范数的函数完全相同例4.9 求矩阵A的三种范数。命令如下： A=[17,0,1,0,15;23,5,7,14,16;4,0,13,0,22;10 ,12,19,21,3;11,18,25,2,19]; a1=norm(A,1) %求A的1—范数 a2=norm(A) %求A的2—范数 ainf=norm(A,inf) %求A的∞—范数
2 2
基于龙格－库塔法，MATLAB提供了求常微分方程数值解的函数，一般调用格式为： [X,Y]=ode45(f,[x0,xn],y0) 其中，X、Y是两个向量， X对应自变量x在求解区间[x1，xn]的一组采样点，其采样密度是自适应的，无需指定； Y是与X对应的一组解， f是一个函数， [x0，xn]代表自变量的求解区间，y0=y(x0)，由方程的初值给定。函数在求解区间[x0，xn]内，自动设立采样点向量X，并求出解函数y在采样点X处的样本值。
4. 多项式的导函数对多项式求导数的函数是： p=polyder(P) 求多项式P的导函数 p=polyder(P,Q) 求P*Q的导函数 [p,q]=polyder(P,Q) 求P/Q的导函数，导函数的分子存入p，分母存入q。
5 矩阵的条件数和迹
1）矩阵的条件数 MATLAB中，计算矩阵A的3种条件数的函数是： (1)cond(A,1) 计算A的1—范数下的条件数 (2)cond(A)或cond(A,2) 计算A的2—范数数下的条件数 (3)cond(A,inf) 计算A的 ∞—范数下的条件数例4.10 求矩阵X的三种条件数。命令如下： A=[2,2,3;4,5,-6;7,8,9]; C1=cond(A,1) C2=cond(A) C3=cond(A,inf)
3 矩阵的秩
MATLAB中，求矩阵秩的函数是rank(A)。例如，求例4.7中方程组系数矩阵D的秩，命令是： D=[2,2,-1,1;4,3,-1,2;8,5,-3,4;3,3,-2,2]; r=rank(D) r= 4 说明D是一个满秩矩阵。
4 向量和矩阵的范数
1）计算向量3种常用范数的函数 (1)norm(V)或norm(V,2) 计算向量V的2—范数 (2)norm(V,1) 计算向量V的1—范数 (3)norm(V,inf) 计算向量V的∞—范数例4.8 已知V，求V的3种范数。命令如下： V=[-1,1/2,1]; v1=norm(V,1) %求V的1—范数 v2=norm(V) %求V的2—范数 vinf=norm(V,inf) %求∞—范数
2．非线性方程组求解
函数fsolve调用格式为： [X,fval]=fsolve(F,X0) 例4.48 求方程组在(1，1，1)附近的解并对结果进行验证。首先建立方程的函数文件fxyz1.m： function F=F(X) x=X(1);y=X(2);z=X(3); F(1)=sin(x)+y+z^2*exp(x); F(2)=x+y*z; F(3)=x*y*z; 在MATLAB命令窗口，输入命令： X=fsolve('fxyz1',[1,1,1]) %求解X的三个分量x、y、z Y=fxyz1(X) %检验所求结果X是否满足原方程组 norm(Y) %求Y向量的模
例4.1 求向量sin(X)的1~3阶差分。设X由[0，2π] 间均匀分布的10个点组成。命令如下： X=linspace(0,2*pi,10); Y=sin(X); DY=diff(Y); %计算Y的一阶差分 D2Y=diff(Y,2); %计算Y的二阶差分，也可用命令diff(DY)计算 D3Y=diff(Y,3); %计算Y的三阶差分，也可用 diff(D2Y)或diff(DY,2)
1．单变量非线性方程求解 MATLAB中，提供了求解单变量方程的函数 [x,favl]=fzero(f,x0,tol)，该函数采用迭代法计算函数f(x)的一个零点，迭代初值为x0，当两次迭代结果小于tol时停止迭代过程。tol的缺省值是eps。注意，函数f(x) 为字符型、内敛函数、匿名函数和 M函数文件的函数句柄。
例4.34设有两个多项式，计算： (1)求f(x)+g(x)、f(x)-g(x)。 (2)求f(x)· g(x)、f(x)/g(x)。在MATLAB命令窗口，输入命令： f=[3,-5,2,-7,5,6];g=[3,5,-3];g1=[0,0,0,g]; f+g1 %求f(x)+g(x) f-g1 %求f(x)-g(x) conv(f,g) %求f(x)*g(x) [Q,r]=deconv(f,g) %求f(x)/g(x)，商式送Q，余式送r。
第四章
4.1 4.2 4.3 4.4
数值计算
数值微积分矩阵和代数方程概率分布和统计分析多项式计算
4.1.1 数值微分 MATLAB中，没有直接提供求数值导数的函数，只有计算向前差分的函数。 DX=diff(X) 计算向量X的向前差分，DX(i)=X(i+1)X(i)，0<i<n。 DX=diff(X,n) 计算X的n阶向前差分， diff(X,2)=diff(diff(X))。 DX=diff(A,n,dim) 计算矩阵A的n阶差分，dim=1时 (缺省状态)，按列计算差分，dim=2，按行计算差分。
4.1.2数值积分 (1)被积函数是一个解析式函数quad(f,a,b,tol)用于求被积函数f(x)在[a,b] 上的定积分，tol是计算精度，缺省值是0.001。注意，调用quad函数时，先要建立一个描述被积函数 f(x)的函数文件或语句函数;f(x)可以是字符型、内敛对象、匿名函数和M函数文件的函数句柄。当被积函数f含有一个以上的变量时， quad函数的调用格式为： quad(f,a,b,tol,g1,g2) 其中 f,a,b,tol等参数的含义同前。
4. 多项式的四则运算 (1)多项式的加减法当两个多项式的次数不同时，要在一个较低次幂的多项式系数向量前补0，使两个系数向量等长。 (2)多项式的乘法函数conv(P1,P2)用于求多项式P1和P2的乘积。 (3)多项式的除法函数[Q,r]=deconv(P1,P2)用于对多项式P1和P2作除法运算。其中Q返回多项式P1除以P2的商式，r返回P1 除以P2的余式。这里，Q和r仍是多项式系数向量。 deconv是conv的逆函数，即有P1=conv(P2,Q)+r。
(3)二重积分
MATLAB提供了计算二重积分的函数： dblquad(f,a,b,c,d,tol) 该函数求f(x,y)在[a,b]×[c,d]区域上的二重积分。参数tol的用法与函数quad完全相同。 1 1 e x y dxdy 例求积分 0 0 命令如下： g=inline(‘exp(-x.^2-y.^2)’); %构造内联函数 dblquad(g,0,1,0,1) %直接调用二重积分函数求解
例4.47 求f(x)=x-+5 在x0=-5和x0=1作为迭代初值时的零点。先编制一个函数文件fz.m： function f=f(x) f=x-1/x+5; 然后，在MATLAB命令窗口，输入命令： fzero('fz',-5) %以-5作为迭代初值 >>Zero found in the interval: [-4.8, 4.2]. fzero('fz',1) >>ans = 0.1926
例4.11 用2种不同的格式求A的特征值和特征向量。命令如下： A=[1,2,2;1,-1,1;4,-12,1]; E=eig(A) [V,D]=eig(A)
4.2.2 线性方程组求解
1．利用左除运算符的直接解法对于线性方程组Ax=b，可以利用左除运算符“\”求解： x=A\b
例4.13 用直接解法求解下列线性方程组。命令如下： A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,1,-4]; b=[13,-9,6,0]'; x=A\b X=inv(A)*b
2 方阵的行列式求方阵A所对应的行列式的值的函数是det(A)。例4.7 用克莱姆(Cramer)方法求解线性方程组。程序如下： D=[2,2,-1,1;4,3,-1,2;8,5,-3,4;3,3,-2,2]; %定义系数矩阵 b=[4;6;12;6]; %定义常数项向量 D1=[b,D(:,2:4)]; %用方程组的右端向量置换 D的第1列 D2=[D(:,1:1),b,D(:,3:4)]; %用方程组的右端向量置换D的第2列 D3=[D(:,1:2),b,D(:,4:4)]; %用方程组的右端向量置换D的第3列 D4=[D(:,1:3),b]; %用方程组的右端向量置换 D的第4列 DD=det(D);x1=det(D1)/DD;x2=det(D2)/DD;x3=det(D 3)/DD; x4=det(D4)/DD;[x1,x2,x3,x4]