Matlab 科学计算

18
>>x1=[0.5 :0.05 :3.0]; >>y1=a(1)*x1.^2+a(2)*x1+a(3); >>plot(x,y,’*’) >>hold on >>plot(x1,y1,’-r’)
19
计算方法中的求积公式
• 微积分是理工科大学生踏入大学校门后首 先要接触到的一门重要基础课程。 • 实践证明，微积分运算也是实际工作中广 泛应用的基本工具，它是科学计算最重要 的组成部分
பைடு நூலகம்
26
trapz实例
• Z = TRAPZ(X,Y,DIM) or TRAPZ(Y,DIM) integrates across dimension DIM of Y. • Example: If Y = [0 1 2 • 3 4 5] • then trapz(Y,1) is [1.5 2.5 3.5] and trapz(Y,2) is [2 8];
10
分段线性插值法
• 高次插值多项式可能产生极大的误差（这 一现象称为荣格（Runge）现象）因此不 宜用太多的点来做插值多项式。 • 为了解决荣格现象，引入了分段线性插值， 即通过n+1个插值点用折线段连接起来逼 近原曲线，这也是计算机绘制图形的基本 原理。
11
荣格（Runge）现象
x=[-5:1:5]; y=1./(1+x.^2); x0=[-5:0.1:5]; y0=lagrange(x,y,x0); y1=1./(1+x0.^2); plot(x0,y0,'--r'); hold on plot(x0,y1,'-b');
9
Lagrange插值法
• Lagrange插值公式如下：
n x xj y ( x) y k k 0 0 xk x j jj k
n

• 对任意k(k=0,…,n), Lagrange插值基函数：
( x x0 )(x x1 )( x xk 1 )(x xk 1 )( x xn ) l k ( x) ( xk x0 )(xk x1 )( xk xk 1 )(xk xk 1 )( xk xn )
17
Matlab拟合函数
• 利用polyfit函数进行多项式拟合 • 格式：p=polyfit(X,Y,N)
p(1)*X^N+p(2)*X^(N-1)+…+p(N)*X+p(N+1)
• 实例： >>x=[0.5 1.0 1.5 2.0 2.5 3.0]; >>y=[1.75 2.45 3.81 4.80 8.00 8.60]; >>a=polyfit(x,y,2) a= 0.5614 0.8287 1.1560
27
>>x=[0:2*pi/100:2*pi]; >>y=sin(x); >>z=trapz(y) z= 0 >> x=[0:pi/100:pi]; >> y=sin(x); >>z=trapz(y)*pi/100 z= 1.9998
trapz实例
28
思考。。。
• 如何利用trapz求下列积分的值
13
插值实例
x = 0:10; y = sin(x); xi = 0:.25:10; yi = interp1(x,y,xi); plot(x,y,'o',xi,yi)
14
插值实例
x=[-5:1:5]; y=1./(1+x.^2); x0=[-5:0.1:5]; y0=interp1(x,y,x0); plot(x0,y0,'*m') hold on plot(x,y,‘bo')
•
25
MATLAB求积函数trapz
• trapz: trapezoidal numerical integration(梯形数值 积分) • Form：Z = TRAPZ(Y) computes an approximation of the integral of Y via the trapezoidal method (with unit spacing). To compute the integral for spacing different from one, multiply Z by the spacing increment. • Comments:For vectors, TRAPZ(Y) is the integral of Y. For matrices, TRAPZ(Y) is a row vector with the integral over each column.
16
拟合
• 在实际中，往往需要从一组实验数据(xi,yi), 中找到自变量和因变量的函数关系y=f(x)。 由于观测数据往往不够准确，因此不要求 y=f(x)经过所有的点(xi,yi),而只要求在给 定点xi上的误差det(i)=f(xi)-yi按某种要求 达到最小。 • 最小二乘拟合采用的策略是使误差的平方 和达到最小。
4
插值和拟合示意图
插值示意图
拟合示意图
5
插值与拟合的应用背景
• 在生产实际和科学研究中，会遇到大量函数，其 中相当一部分是通过测量和实验观测得到的，虽 然其函数关系y=f(x)在某区间上客观存在，但却 不知道具体的解析表达式。只能得到函数在区间 [a,b]上的一些离散点的函数值、导数值等。 • 我们希望对这样的函数用一个比较简单的函数的 表达式来近似地给出其整体描述。有时是虽然有 明确表达式，但函数本身过于复杂而不便于进行 数值计算，同样希望构造一个既能反映函数性质 又便于计算的简单函数，来近似代替原来的函数， 插值就是寻求近似函数的方法之一。
3
e
0
0.5t
sin(t / 6)dt
29
解线性方程组的应用背景
• 求解线性方程组在科学计算中最常遇到， 数值方法如三次样条、有限元法等都要涉 及到求解线性方程组。 • 求解线性方程组的实际应用也是非常广泛： 天气预报、石油勘探、地质测绘、电路分 析、分子结构预测等。
30
解线性方程组的理论方法
Matlab与科学计算
计算机学院 刘咏梅 Email:liuyongmei@
1
第二章 MATLAB数值计算基础
• Introduction to MATLAB • MATLAB Basics of Numerical Computing • MATLAB Programming(branching and Loops) • MATLAB Programming(Predefined Functions) • MATLAB Graphic and Image Processing • Scientific Computing using MATLAB
24
复合梯形求积公式
• 复合梯形求积公式的由来：为了提高精度，常 把区间[a,b]分成若干个小区间，在每个小区间 上使用次数不高的Newton-Cotes求积公式， 如梯形公式和Simpson公式，再将所有的子区 间上的积分值累加起来。 复合梯形求积公式：f(x)在区间[a,b]上的积分 值为
n 1 ba [ f (a ) 2 f ( xk ) f (b)] 2n k 1
b a k 0 n
Cotes系数
22
Cotes系数
当n=1时，c0=1/2, c1=1/2
梯形公式
Simpson公式
当n=2时，c0=1/6, c1=2/3, c2=1/6 ………
当n=3时，c0=1/8, c1=3/8, c2=3/8, c3=1/8
23
梯形公式
In
b a
ba f x dx ( f (b) f (a)) 2
20
计算方法中的求积公式
• Newton-Cotes系列数值求积公式 梯形公式 复合梯形公式
21
Newton-Cotes系列数值求积公式
• 一种等间隔闭式求积公式 • 将积分区间[a,b]划为n等分，步长 h=(b-a)/n，选取等距离节点xk=a+kh，构造下面 的插值型求积公式：
I n f x dx (b a) C k( n ) f xk
• 克莱姆法则（Cramer’s Rule）：若线性方程组⑴的系数 行列式 D≠0，则线性方程组⑴有唯一解，其解为 xj=Dj/D，其中Dj是把D中第j列元素对应地换成常数项， j=1,…,n,而其余各列保持不变所得到的行列式。
31
解线性方程组的理论方法
• 如果按克莱姆法则求解，假设计算机每秒可计算1亿次乘 法，共n+1个行列式，每个行列式计算需要n!(n+1)次乘 法，故共需(n2-1) n!次乘法。当问题规模n=10时的计算 量耗时3.59秒，而n=100时已经非常大了，不能接受。 • 因此说，理论上比较完善的方法在实际中却用处不大。
2
本章学习的目标：
• 利用MATLAB进行插值和拟合 • 利用MATLAB求积分和微分 • 利用MATLAB求解线性方程组
3
插值和拟合
• 插值：利用函数f(x)在某区间中若干个点的 函数值，做适当特定的函数f’(x)，使区间 上其它点的值用f’(x)作为f(x)的近似。 • 拟合：已知某函数f(x)若干离散值，通过调 整f(x)中若干待定系数，使f(x)与已知点集 的差别最小。
15
二维插值函数interp2
• ZI = interp2(X,Y,Z,XI,YI,’method’) interpolates to find ZI. • The values of the underlying 2-D function Z at the points in matrices XI and YI. • Matrices X and Y specify the points at which the data Z is given. • Method=nearest/linear(default)/spline/cub ic