实验四 数值运算

3. 数值微分
(1) 差分使用diff 函数实现。 >>x=1:2:9
>>diff(x)
(2) 可以用因变量和自变量差分的 结果相除得到数值微分。 >>x=linspace(0,2*pi,100); >>y=sin(x);
>>plot(x,y)
>>y1=diff(y)./diff(x); >>plot(x(1:end-1),y1)
2. 多项式插值和拟合 有一组实验数据如下表所示。
x 1 2 3 4 5 6 7 8 9 10
y
16
32
70
142
260
436
682
1010 1432 1960
请分别用拟合(二阶至三阶)和插值(线性和三次样条)的方法来估 测 X=9.5 时Y 的值。 以下是实现一阶拟合的语句。
>>x=1:10 >>y=[16 32 70 142 260 436 682 1010 1432 1960] >>p1=ployfit(x,y,1) %一阶拟合 >>y1=ployval(p1,9.5) %计算多项式p1 在x＝9.5 时的值
ห้องสมุดไป่ตู้
例题
x=0:0.1:1; y = [-0.447, 1.978,3.28, 6.16, 7.08, 7.34, 7.66, 9.56, 9.48, 9.30, 11.2] a1 = polyfit(x, y, 1); x1 = linspace(0,1); y1 = polyval(a1,x1) figure; subplot(2,2,1); plot(x,y,'o',x1,y1,'b'); a2 = polyfit(x,y,2); y2 = polyval(a2,x1); subplot(2,2,2); plot(x,y,'o',x1,y2,'m'); a3 = polyfit(x,y,3); y3 = polyval(a3,x1); subplot(2,2,3); plot(x,y,'o',x1,y3,'r'); a9 = polyfit(x,y,9); y9 = polyval(a9,x1); subplot(2,2,4); plot(x,y,'o',x1,y9,'c');
4. 数值积分
(1) cumsum 函数求累计积分， (3). 辛普森数值积分 q=quad('f',a,b) 表示使用自适应递归 >>x=ones(1,10) 的辛普森方法从积分区间a 到b 对函 数f(x)进行积分，积分的相对误差在 >> cumsum(x) 1e-3 范围内。输入参数中的'f'是一个 (2)trapz 函数用梯形法求定积分，字符串，表示积分函数的名字。当输 即曲线的面积。 入的是向量时，返回值也必须是向量 形式。 >> x=linspace(0, pi,100); >> y=sin(x); >> S=trapz(y,x)
4，利用梯形法和辛普森法求定积分 的值， 并对结果进行比较。如果积分区间改为-5～5 结果有 何不同？梯形积分中改变自变量x 的维数，结果有 何不同？
练习5：图A1 是瑞士地图，为了算出其国土面积，首先对地图作如下测量：
以由西向东方向为X 轴，由南到北方向为Y 轴，选择方便的原点，并将从最西 边界点到最东边界点在X 轴上的区间适当划分为若干段，在每个分点的Y 方向 测出南边界点和北边界点的Y 坐标Y1 和Y2，这样就得到了表1，根据地图比例 尺知道18mm 相当于40km，试由测量数据计算瑞士国土近似面积，与其精确值 41228km2 比较。地图的数据见下表1-2。
• 格式：zi=interp2(x,y,z,xi,yi,’method’)
x,y为已知的点坐标向量，z为矩阵（x,y对应点的值）xi，
yi 为插入点的 X,Y 坐标向量，‘method’：同上
%立体插值 x=(-4:1:4);y=x; [x1,y1]=meshgrid(x,y); z=peaks(x1,y1); subplot(2,1,1); mesh(x1,y1,z); xi=(-4:0.2:4); yi=xi'; zi=interp2(x,y,z,xi,yi,'cubic'); subplot(2,1,2); mesh(xi,yi,zi+20);
练习
1. 用函数roots 求方程x2 − x −1 = 0的根。 2. y = sin x,0≤x≤2π，在n 个节点(n 不要太大，如取 5～11)上用线性和三次样条插值方法，计算m 个插 值点(m 可取50～100)的函数值。通过数值和图形输 出，将两种插值结果与精度值进行比较。适当增加n， 再作比较。
练习
3. 大气压强p 随高度x 变化的理论公式为p =1.0332e−( x+500) / 7756，为验证这一公式，测得某地大气 压强随高度变化的一组数据如表所示，试用插值法和 拟合法进行计算并绘图，看哪种方法较为合理，且总 误差最小
高度 /m 压强 /Pa 0 0.9689 300 0.9322 600 0.8969 1000 0.8519 1500 0.7989 2000 0.7491
实验四 数值运算
一、实验目的 掌握 MATLAB 的数值运算及其运算中所用到的函数。 二、实验内容： (1) 多项式运算。 (2) 多项式插值和拟合。 (3) 数值微积分。
(3) 多项式的乘、除法分别用函数conv 和deconv 实现
>>S1=[ 2 3 11 ] >>S2=[1 3 -5 4 7 ] >>S3=conv(S1,S2) >>S4=deconv(S3,S1)
提示：由高等数学的知识可知，一条曲线的 定积分是它与x 轴所围成的面积，那么两条 曲线所围成的面积可由两条曲线的定积分相 减得到。
例4-11 某观测站测得某日6:00时至18:00时之间每隔2小时的 室内外温度(℃)，用3次样条插值分别求得该日室内外6:30 至17:30时之间每隔2小时各点的近似温度(℃)。 设时间变量h为一行向量，温度变量t为一个两列矩阵，其中 第一列存放室内温度，第二列储存室外温度。命令如下： h =6:2:18; t=[18,20,22,25,30,28,24;15,19,24,28,34,32,30]'; XI =6.5:2:17.5 YI=interp1(h,t,XI,‘spline’) %用3次样条插值计算
例： 假设有一个汽车发动机在转速为 2000r/min 时， 温度(单位为℃)与时间(单位为s)的5 个测量值如表所示。
>> t=[0 1 2 3 4 5]'; % 输入时间 y=[0 20 60 68 77 110]'; % 输入温度 y1=interp1(t,y,2.5) % 要内插的数据点为 2.5 y1=interp1(x,y,[2.5 4.3]) %内插数据点为2.5,4.3，注意采用[ ]放入多个内插点 y1=interp1(x,y,2.5,'cubic') %以三次方程式对数据点2.5 作内插 y1=interp1(x,y,2.5,'spline') %以spline 函数对数据点 2.5 作内插
(4) 多项式求根用函数roots
>> S1=[ 2 4 2 ] >> roots(S1)
(5) 多项式求值用函数polyval
>> S1=[ 2 4 1 -3 ] >> polyval(S1,3) %计算x＝3 时多项式的值 >> x＝1:10 >> y=ployval(S1,x) %计算x 向量对应的值得到y 向量
4. 数值积分
例 计算二重定积分 (1) 建立一个函数文件fxy.m： function f=fxy(x,y) global ki; ki=ki+1; %ki用于统计被积函数的调用次数 f=exp(-x.^2/2).*sin(x.^2+y); (2) 调用dblquad函数求解。 global ki;ki=0; I=dblquad('fxy',-2,2,-1,1) ki I= 1.57449318974494 ki = 1038