实验数学五：MATLAB在数值计算中的应用

mesh(x0,y0,temps) 结果见图5.6
分别用线性性插值和三次样条插值求已知点的温度。
输入命令 ：
19
>> t=interp2(x0,y0,temps,[1.5 2 2.5 3.5],[1.5 1.6 2 4.5],'liner')
t= 76.2500 70.2000 62.0000 NaN
>> T=interp2(x0,y0,temps,[1.5 2 2.5 3.5],[1.5 1.6 2 4.5],'spline')
T= 71.4531 65.5200 60.9688 188.8906
（2）以平滑数据,在x、y方向上每隔0.2个单位的地方 进行插值画出线性和三次样条插值的温度分布曲面图.
输入以下命令得温度的线性插值曲面图: >> x=1:0.2:5; >> y=1:0.2:3; >> z=interp2(x0,y0,temps,x',y,'linear'); >> mesh(x,y,z)
20
>> xlabel('x轴'); >> ylabel('y轴'); >> zlabel('z轴') >> title('线性插值曲面图') 结果见图5.7
而method为可选参数，有四种选择： ‘nearest’ 最邻近插值 ‘linear’ 线性插值
‘spline’ 三次样条插值
‘cubic’ 三次插值
缺省时, 是线性插值 17
例8：测得平板表面3×5网格点处的温度分别为： 82 81 80 82 84 79 63 61 65 81 84 84 82 85 86
且x1 x2 xn ; y1 y2 yn ,
求点(x*, y*)( (xi , yj ))处的插值z*.
对上述问题，MATLAB提供了二维插值函数 interp2,其基本格式为：
16
z=interp2(x0,y0,z0,x,y,’method’)
其中x0,y0是自变量。X0,y0的分量值必须是单调 递增的。X0和y0分别是m维和n维向量，分别表示已 知数据点的横、纵坐标向量，z0是m*n维矩阵，标明 相应于所给数据网格点的函数值。向量x,y是待求函数 值所给定网格点的的横、纵坐标向量，x,y的值分别不 能超出x0,y0的范围。
试求在平板表面坐标为（1.5,1.5）,（2,1.6）,（2.5,2） （3.5,4.5）处的温度，并作平板表面的温度分布曲面 z=f(x,y)的图形，
（1）先在三维坐标画出原始数据，画出粗糙的温 度分布曲图.
输入以下命令：
x0=1:5;
y0=1:3;
temps=[82 81 80 82 84;79 63 61 65 81;84 84 82 85 86]; 18
直线，即x轴。
2
方程在[0,10]区间从图中可看出有4个解，分别在 0，3，6，9附近，
所以用命令：
>> f=inline('5*x.^2.*sin(x)-exp(-x)');
>> fsolve(f,[0,3,6,9])
ans =
0.5018 3.1407 6.2832 9.4248
3
2、零点法 格式：fzero(‘f’,x0): 求函数f在x0 附近的零点。
>> t=interp1(hours,temps,[3.2 6.5 7.1 11.7])
%线性插值
t=
10.2000 30.0000 30.9000 24.9000
13
>> T=interp1(hours,temps,[3.2 6.5 7.1 11.7],'spline') %三次样条插值
T= 9.6734 30.0427 31.1755 25.3820
9
由图可看出可用二次多项式拟合。
再输入命令 ：
>> p=polyfit(x,y,2)
p= 0.5614 0.8287 1.1560
即二次拟合多项式为
f (x) 0.5614x2 0.8287x 1.1560
10
画出离散点及拟合曲线： 输入命令 ： >> x1=0.5:0.05:3.0; >> y1=polyval(p,x1); >> plot(x,y,'*r',x1,y1,'-b') 结果见图5.4
最常用的ixiy曲线拟合是最小二乘法曲线拟合曲线拟合是最小二乘法曲线拟合拟合结果可使误差的拟合结果可使误差的上尽可能8平方和最小即找出使21fxniiiy??最小的fx
MATLAB在数值计算中的应用
5.1实验目的
在工程技术中，大量的实际问题都需要进行近似 处理，从而产生不同问题的数值计算方法。而 MATLAB具有强大的数值运算功能，本实验的目的是 学会用MATLAB软件进行一些数值运算，包括代数方 程求根、插值问题和曲线拟合问题等。
>> x1=fsolve(f,0)
x1 = 0.5671
例2 求方程5x2sinx-e-x =0的根。
先画图观察根的个数及大概位置。
输入命令 ： >> fplot('[5*x^2*sin(x)-exp(-x),0]',[0,10])
结果见图5.1
注意，[5*x^2*sin(x)-exp(-x),0]中的 […,0]是作y=0
>> ylabel('温度')
14
>> gtext('线性插值') >> gtext('三次样条插值')
结果见图5.5
15
四、二维插值 对二维插值问题，MATLAB分别给出了针对插值基
点为网格节点的插值函数及针对插值基点为散乱节点的 插值函数调用格式。
1、 用MATLAB作网格节点数据的插值 已知m n个节点：(xi,yj,zij),(i=1,2, m; j 1, 2, , n)
22
例 9 山区地貌：
在某山区测得一些地点的高程如下表5.2。平面区 域为 1200<=x<=4000,1200<=y<=3600) 试作出该山区的地貌图和等高线图，并对几种插值方 法进行比较。
Yx
0 400 800 1200 1600 2000 2400 2800 3200 3600 4000 4400 4800
x=
[ 6.7]
[ -1.7]
7
>> y=vpa(y,2)
y= [ .14] [ -.60]
二、曲线拟合
已知离散点上的数据集 [(x1, y1)(x2 , y2 ) (xn , yn )],
求得一解析函数y=f(x)使y=f(x)在原离散点 xi 上尽可能
接近给定 yi 的值，这一过程叫曲线拟合。最常用的
21
再输入以下命令得温度的三次样条插值曲面图: >> z=interp2(x0,y0,temps,x',y,'spline'); >> mesh(x,y,z) >> xlabel('x轴'); >> ylabel('y轴'); >> zlabel('z轴') >> title('三次样条插值曲面图') 结果见图5.7
5.2实验内容
一、代数方程求根
代数方程求根有各种近似处理方法，下面给出
MATLAB两种常用的调用格式： 最小二乘法
格式：fsolve(‘f’,x0):求方程f=0在估计值x0 附近
的近似解。
1
例1 求方程x-e-x 0在x 0附近的根。 解 输入命令 ： >> f=inline('x-exp(-x)');
曲线拟合是最小二乘法曲线拟合，拟合结果可使误差的
平方和最小，即找出使
n
f (xi ) yi 2
i 1
最小的f(x).
8
格式：p=polyfit(x,y,n).
说明：求出已知数据x,y 的n次拟合多项式f(x)的系 数p,x 必须是单调的。
例6 已知某函数的离散值如表5.1
xi
0.5
1.0
0 400 800 1200 1600 2000 2400 2800 3200 3600 4000 4400 4800 5200 5600
370 470 550 600 670 690 670 620 580 450 400 510 620 730 800 850 870 850 780 720 650 500 650 760 880 970 1020 1050 1020 830 900 700 300 740 880 1080 1130 1250 1280 1230 1040 900 500 700 830 980 1180 1320 1450 1420 1400 1300 700 900 850 880 1060 1230 1390 1500 1500 1400 900 1100 1060 950 910 1090 1270 1500 1200 1100 1350 1450 1200 1150 1010 950 1190 1370 1500 1200 1100 1550 1600 1550 1380 1070 1430 1430 1460 1500 1550 1600 1550 1600 1600 1600 1550 1420 1430 1450 1480 1500 1550 1510 1430 1300 1200 980 1380 1410 1430 1450 1470 1320 1280 1200 1080 940 780 1370 1390 1410 1430 1440 1140 1110 1050 950 820 690 1350 1370 1390 1400 1410 960 940 880 800 690 570
例7 在1-12的11小时内，每隔1小时测量一次温度，测 得的温度依次为：5，8，9，15，25，29，31，30， 22，25，27，24。试估计在 3.2,6.5,7.1,11.7小时的温 度值。
解 输入命令 ：
>> hours=1:12; >> temps=[5 8 9 15 25 29 31 30 22 25 27 24];
>> [x,y]solve('x*y=1','x-11*y=5')
x=
[ 5/2+1/2*69^(1/2)]
[ 5/2-1/2*69^(1/2)]
y=
[ -5/22+1/22*69^(1/2)]
[ -5/22-1/22*69^(1/2)]
如果化成数值解，用命令vpa 如上例：
>> x= vpa(x,2)
fzero(‘f’,[x1,x2]): 求函数f在区间[x1,x2] 上唯一零点。
例3 求方程x2-4x-5=0的根。
先画图观察根的个数及大概位置。 输入命令 ： >> fplot('[x^2-4*x-5,0]',[-10,10]) 结果见图5.2
4
从图中可看出方程在[-2,0]及[4，6]区间上各有一根，
11
三、一维插值
已知离散点上的数据集[(x1, y1)(x2 , y2 ) (xn , yn )],
求得一解析函数连接自变量相邻的两个点，并求得两点
间的数值，这一过程叫插值。
MATLAB在一维插值函数interp1中，提供了四种 插值方法选择：线性插值、三次样条插值、立方插值 和最近邻点插值。interp1的本格式为：
比较发现，两种结果有差异，这是因为插值是一 个估计或猜测的过程。
两种插值的画图如下；
输入命令 ：
>> t0=1:0.1:12;
>> T0=interp1(hours,temps,t0,'spline');
>> plot(hours,temps,'+',t0,T0,hours,temps,'r:')
>> xlabel('时间');
例4 求方程ax2 bx c 0的解
解 输入命令 ： >> solve('a*x^2+b*x+c') ans =
[ 1/2/a*(-b+(b^2-4*a*c)^(1/2))] [ 1/2/a*(-b-(b^2-4*a*c)^(1/2))]
6
例5
求方程组
xy=1 x 11y

5
解 输入命令 ：
yi=interp1(x，y，xi，'method') 其中x,y分别表示数据点的横、纵坐标向量，x 必须
单调，xi为需要插值的横坐标数据（或数组），xi不能 超出x的范围，而method为可选参数，有四种选择：
‘nearest’ ：最邻近插值 ‘linear’ ： 线性插值；
12
‘spline’ ： 三次样条插值； ‘cubic’ ： 立方插值。 缺省时： 分段线性插值。
再输入命令 ：
>> x1=fzero('x^2-4*x-5',[-2,0]) x1 = -1
>> x2=fzero('x^2-4*x-5',[4,6])
x2 =
5
5
3、代数方程的符号解 格式：
solve('f',):求代数方程f=0的根； solve('eqn1','eqn2',...,'eqnN'):求n个代数方程的根；
1.5
2.0
2.5
3.0
yi
1.75
2.45
3.81
4.80
7.00
8.65
求二次拟合多项式.
先画函数离散点的图形
输入命令 ：
>> x=[0.5 1.0 1.5 2.0 2.5 3.0];
>> y=[1.75 2.45 3.81 4.80 7.00 8.60];
>> scatter(x,y,5)
结果见图5.3