MATLAB软件使用简介.doc

Matlab软件使用简介
一、基本内容
1. 向量的产生
基本格式：start: step: end（step缺省时为1）
x=1:5
x=0:0.1:2*pi
dot 向量点积 cross 向量叉积 .* 对应元素乘积
Matlab的每条命令后，若为逗号或无标点符号，则显示命令的结果；若命令后为分号，则禁止显示结果。
linspace 线性等分向量 a=linspace(0,2*pi,100)
logspace 对数等分向量 a=logspace(0,2*pi,100)
100为插入点数（99等分），默认值为100.
2. 矩阵的输入
约定：（1）元素之间用空格或逗号隔开；
（2）用中括号（[ ]）把所有元素括起来；
（3）用分号（；）说明行结束。
矩阵输入时，按Enter键表示开始输入新的一行，输入矩阵时，要求所有的行具有相同的列。
例如：a=[1 2 3; 4,5 6; 7 8,9]
该矩阵一直保存在工作空间，直至被修改。
x=1:6
y=sin(x)
z=cos(x)
b=[x;y;z]
特殊矩阵：
a=[ ] 产生一个空矩阵
b=zeros(m,n) 产生一个m行、n列的零矩阵
c=ones(m,n) 产生一个m行n列的元素全为1的矩阵
d=eye(m,n) 产生一个m行、n列的单位矩阵

3. 大矩阵中抽取小矩阵
c1=b(:,1:2); c2=b(:,5:6) ; d=[c1,c2]
d=[c1;c2]
或d=b(:,[1:2,5:6])

4. 固定变量
(1) ans: 在没有定义变量名时，系统默认变量名为ans；
(2) eps: 容许误差，非常小的数；
(3) pi: 即圆周率；
(4) i: 虚数单位；
(5) inf: 表示正无穷大，由1/0运算产生；
(6) NaN: 表示不定值，由inf/inf或0/0运算产生。

5. 基本运算
5.1 算算术运算符
+ 加 - 减 * 矩阵乘法 .* 数组乘（对应元素相乘）
^ 矩阵幂 .^ 数组幂（各个元素求幂）
./ 数组除（对应元素除）
\ 左除或反斜杠 / 右除或斜杠
如果a为一个非奇异矩阵：
a\b = inv(a)*b： 表示a*x=b的解；
a/b = b*inv(a)： 表示项x*a=b的解。
例如：a=[1 2 3; 4 2 6; 7 4 9]
b=[4;1;2]
a\b
5.2 2.2 关系运算符
= = 等号 ~= 不等号 < 小于 > 大于
<= 小于或等于 >= 大于或等于

5.3 逻辑运算符
& 逻辑与 | 逻辑或 ~ 逻辑非 xor 异或
any 有非零元则为真
all 所有元素均非零则为真
6. 矩阵的基本操作
' 转置”
inv 矩阵求逆
det 行列式的值
[v d]=eig(a) 特征值和特征向量
rank 秩
trace 迹
size 矩阵的行数和列数
diag 对角矩阵和矩阵对角线
fliplr 从左自右翻转矩阵
flipud 从上到下翻转矩阵
roy90 矩阵翻转90度
tril 矩阵的下三角
triu 矩阵的上三角
7. 常用函数
（1）clc 清除指令窗口
（2）clear 从内存中清除变量和函数
（3）who 列出工作内存中的变量名

（4）whos 列出工作内存中的变量详细信息
（4）pause 暂停
二、绘图函数
1. plot是绘制一维曲线的基本函数，但在使用此函数之前，我们需先定义曲线上每一点的x 及y坐标。
下例可画出一条正弦曲线：
x=linspace(0, 2*pi, 100); % 100个点的x坐标
y=sin(x); % 对应的y坐标
plot(x,y);
若要画出多条曲线，只需将坐标对依次放入plot函数即可
plot(x, sin(x), x, cos(x));
若要改变颜色，在坐标对后面加上相关字串即可：
plot(x, sin(x), 'c', x, cos(x), 'g');
若要同时改变颜色及图线型态（Line style），也是在坐标对后面加上相关字串即可：
plot(x, sin(x), 'co', x, cos(x), 'g*')
plot绘图函数的参数
符号
颜色
符号
线形
符号
线形
y
黄色
.
点
--
虚线
m
洋红色
o
圆
d
菱形
c
青色
×
叉号
>
向右三角形
r
红色
+
加号
<
向左三角形
g
绿色
*
星号
s
正方形
b
蓝色
-
实线
P
正五角星
w
白色
:
点线

h
正六角星
k
黑色
-.
点划线

图形完成后，我们可用axis([xmin,xmax,ymin,ymax])函数来调整图轴的范围： axis([0, 6, -1.2, 1.2]);
此外，Matlab也可对图形加上各种注解与处理：
xlabel('Input Value'); % x轴注解
ylabel('Function Value'); % y轴注解
title('Two Trigonometric Functions'); % 图形标题
legend('y = sin(x)','y = cos(x)'); % 图形注解
grid on; % 显示格线
注意：如果标题太长，需要分行，请用如下格式
str={'你需要多少行？';'3行可以吗';'试试就知道了'}
title(str)

我们可用subplot来同时画出数个小图形于同一个视窗之中：
subplot(2,2,1); plot(x, sin(x));title('y=sinx') 在2*2矩阵中的第一个图像
subplot(2,2,2); plot(x, sin(2*x)); title('y=sin2x')
subplot(2,2,3); plot(x, sin(3*x)); title('y=sin3x')
subplot(2,2,4); plot(x, sin(4*x)); title('y=sin4x')

Matlab还有其他各种二维绘图函数，以适合不同的应用。
bar 长条图
errorbar 图形加上误差范围
fplot 较精确的函数图形
polar 极坐标图
hist 直方图
rose 绘制角度直方图（玫瑰图）
stairs 阶梯图
stem 针状图
fill 多边形填充图
feather 羽毛图
compass 罗盘图
quiver 向量场图
以下我们针对每个函数举例。
（1）当数据点数量不多时，长条图是很适合的表示方式：
close all; % 关闭所有的图形视窗
x=1:10;（无特殊说明步长为1） y=rand(size(x)); bar(x,y);

（2）如果已知数据的误差量，就可用errorbar来表示。下例以单位标准差来做资的误差量：
x = linspace(0,2*pi,30);
y = sin(x);
e = std(y)*ones(size(x)); % std求标准差
errorbar(x,y,e)

（3）对於变化剧烈的函数，可用fplot来进行较精确的绘图，会对剧烈变化处进行较密集的取样，如下例：
fplot('sin(1/x)', [0.02 0.2]); % [0.02 0.2]是绘图范围

（4）若要产生极座标图形，可用polar：
theta=linspace(0, 2*pi);
r=cos(4*theta);
polar(theta, r);

（5）对於大量的数据，我们可用hist来显示资料的分布情况和统计特性。下面几个命令可用来验证由randn产生的正态分布：
x=randn(5000, 1); % 产生5000个均值为0，方差为1的正态分布随机数
hist(x,20); % 20代表长条的个数

  （6）rose和hist很接近，只不过是将数据大小视为角度，数据个数视为距离，并用极座标绘制表示：
x=randn(1000, 1);
rose(x);

（7）stairs可画出阶梯图：
x=linspace(0,10,50);
y=sin(x).*exp(-x/3);
stairs(x,y);

（8）stems可产生针状图，常被用来绘制数位讯号：
x=linspace(0,10,50);
y=sin(x).*exp(-x/3);
stem(x,y);

（9）fill将资料点视为多边行顶点，并将此多边行涂上颜色：
x=linspace(0,10,50);
y=sin(x).*exp(-x/3);
fill(x,y,'b'); % 'b'为蓝色

（10）feather将每一个数据点视复数，并以箭号画出：
theta=linspace(0, 2*pi, 20);
z = cos(theta)+i*sin(theta);
feather(z);

（10）compass和feather很接近，但每个箭号的起点都在圆点：
theta=linspace(0, 2*pi, 20);
z = cos(theta)+i*sin(theta);
compass(z);

(11) quiver向量图
2. plot3用来绘制空间曲线
格式：plot3(x,y,z, cs)
其中若x,y,z为长度相等的向量，则在空间表示一条曲线；若x,y,z为m*n的矩阵，则每一列对应一条曲线。cs表示曲线的颜色、线形等性质。
例螺旋线：t=0:pi/50:10*pi; x=sin(t); y=cos(t); plot3(x,y,t)
例多条曲线：x=-3:0.1:3;y=1:0.1:5; [X,Y]=meshgrid(x,y);
Z=(X+Y).^2; plot3(X,Y,Z)

（meshgrid(x,y)的作用是产生一个以向量x为行、向量y为列的矩阵）
3. surf用来绘制空间曲面
格式：surf(x,y,z)
x,y,z为矩阵，显式函数的图形。
例：显示的图形。
n=-8:0.5:8;m=n'; [x y]=meshgrid(n,m);t=sqrt(x.^2+y.^2)+eps;%加eps是因为所选的点中有（0，0）
z=sin(t)./t;surf(x,y,z)
4. mesh用来绘制网格曲面
格式同上
5. scatter离散点绘图
格式：scatter(x,y,cs)
例: x=0:pi/10:2*pi;y=sin(x);scatter(x,y,'+')
6. 其它函数
（1）hold on: 保持当前图形, 以便继续在当前图上画图；
hold off：释放当前窗口
（2）figure(n): 打开或创建图形窗口n
三、微积分内容
1. limit求极限（必须首先声明变量）
格式：limit(F,x,a) 求
limit(F,a) 求
limit(F) 求
limit(F,x,a,'right') 的右极限，即
limit(F,x,a,'left') 的左极限，即
例：syms x; %声明x为符号变量（cosn c
limit((1+1/x)^x,x,inf)；
limit(sin(x)/x,x,0); limit(sin(x)/x,0); limit(sin(x)/x);
syms h; limit((sin(x+h)-sin(x))/h,h,0)
2. diff求导数
格式：diff(F) 求函数的一阶导数
diff(F, v) 求函数对变量一阶偏导数
diff(F,n) 求函数的阶导数
diff(F, v, n) 求函数对变量的阶偏导数
例：diff(x*y*z); diff(x*y*z,2) diff(diff(x*y*z,x),y) 求导对象不同
diff(x^2*cos(x),x,2)
diff('f(x^2)')
3. 积分函数
（1）符号积分
int(f) 传回f对预设独立变数的积分值
int(f,'t') 传回f对独立变数t的积分值
int(f,a,b) 传回f对预设独立变数的积分值，积分区间为[a,b]，a和b为数值式
int(f,'t',a,b) 传回f对独立变数t的积分值，积分区间为[a,b]，a和b为数值式
int(f,'m','n') 传回f对预设变数的积分值，积分区间为[m,n]，m和n为符号式，注意有单引号
例 syms x y int(x) int(x,x) int(x*y,x)
int(x,1,2) int(x,x,1,2)
注意：int(x,x,a,b) int(x,x,'a','b')的区别。
例：计算二重积分,其中由轴，轴及围成。
syms x y
c=int(int(x+y,y,0,1-x),x,0,1)
vpa(c,10) %给出数值型符号结果
（2）数值积分
trapz(x,y) 梯形法沿列方向求函数Y关于自变量X的积分
cumtrapz(x,y) 梯形法沿列方向求函数Y关于自变量X的累计积分
quad(fun,a,b) 采用递推自适应Simpson法计算积分（低阶方法）

quadl(fun,a,b) 采用递推自适应Lobatto法求数值积分（高阶方法）
dblquad(fun,xmin,xmax,ymin,ymax,tol) 矩形区域二重数值积分
triplequad(fun,xmin,xmax,ymin,ymax,zmin,zmax,tol) 长方体区域三重数值积分
例：计算定积分的值。
符号积分：syms x; int(exp(-sin(x)),x,0,4)
数值积分：参看ex01.m
例：计算积分。
quad('sin(x)./x',0,2) quadl('sin(x)./x',0,2)
注意：quad函数和quadl函数求积分时，其中求解变量是以向量来计算的，所以这里的sin(x)和x都是向量，因此要用点除；同时被积函数用' '括起来。
例：计算二重积分，
dblquad('x*y',0,1,0,1)
例：triplequad(' y*sin(x)+z*cos(x)',0,pi,0,1,-1,1)

注意：最新版的2009a有关于一般区域二重积分的计算函数quad2d，如果用的不是2009a，那么你可以利用NIT工具箱里的quad2dggen函数。
quad2dggen(fun,xmin,xmax,ymin,ymax,tol) 任意区域上二重数值积分
如果既没有NIT工具箱用的也不是2009a，可以利用两次quadl函数。

例：计算二重积分，
quadl(@(x) arrayfun(@(xx) quadl(@(y) xx+y,0, xx),x),0,1)
% @创建函数
4. 无穷级数
格式：r = symsum(s) r = symsum(s,v)
r = symsum(s,a,b) r = symsum(s,v,a,b)
syms k n x; symsum(k^2)
symsum(k^2,1,2) symsum(k^2,1,inf)
5. 方程求解
基本函数:
solve(eq)
solveeq,var)
solve(eq1,eq2,...,eqn)
g = solve(eq1,eq2,...,eqn,var1,var2,...,varn)
例 solve('2*x+1=0','x')或syms x; solve('2*x+1=0')
S = solve('x + y = 1','x - 11*y = 5')
S.y
S.x
roots(vec) vec为系数组成的向量
例 roots([2 1])
线性方程的求解分为两类：一类是方程组求唯一解或求特解，另一类是方程组求无穷解即通解。可以通过系数矩阵的秩来判断：
若系数矩阵的秩r=n（n为方程组中未知变量的个数），则有唯一解；
若系数矩阵的秩r>A=[5 6 0 0 0
1 5 6 0 0
0 1 5 6 0
0 0 1 5 6
0 0 0 1 5];
B=[1 0 0 0 1]';
R_A=rank(A) %求秩
X=A\B %求解
运行后结果如下
R_A =
5
X =
2.2662
-1.7218
1.0571
-0.5940
0.3188
这就是方程组的解。
或用函数rref求解：
C=[A,B] %由系数矩阵和常数列构成增广矩阵C
R=rref(C) %将C化成行最简行
R =
1.0000 0 0 0 0 2.2662
0 1.0000 0 0 0 -1.7218
0 0 1.0000 0 0 1.0571
0 0 0 1.0000 0 -0.5940
0 0 0 0 1.0000 0.3188
则R的最后一列元素就是所求之解。
例 求方程组的一个特解。
解：
>>A=[1 1 -3 -1;3 -1 -3 4;1 5 -9 -8];
>>B=[1 4 0]';
>>X=A\B %由于系数矩阵不满秩，该解法可能存在误差。
X =[ 0 0 -0.5333 0.6000]’（一个特解近似值）。
若用rref求解，则比较精确：
>> A=[1 1 -3 -1;3 -1 -3 4;1 5 -9 -8];
B=[1 4 0]';
>> C=[A,B]; %构成增广矩阵
>> R=rref(C)
R =
1.0000 0 -1.5000 07500 1.2500
0 1.0000 -1.5000 -1.7500 -0.2500
0 0 0 0 0
由此得解向量X=[1.2500 – 0.2500 0 0]’（一个特解）。
5.1.2 利用矩阵的LU、QR和cholesky分解求方程组的解
（1）LU分解：
LU分解又称Gauss消去分解，可把任意方阵分解为下三角矩阵的基本变换形式（行交换）和上三角矩阵的乘积。即A=LU，L为下三角阵，U为上三角阵。
则：A*X=b 变成LU*X=b
所以X=U\(L\b) 这样可以大大提高运算速度。
命令 [L，U]=lu (A)
例 求方程组的一个特解。
解：A=[4 2 -1;3 -1 2;11 3 -1];
B=[5 7 14]';
D=det(A)
[L,U]=lu(A)
X=U\(L\B)
见ex04.m
例 求方程组的一个特解。
解：
>>A=[4 2 -1;3 -1 2;11 3 0];
>>B=[2 10 8]';
>>D=det(A)
>>[L,U]=lu(A)
>>X=U\(L\B)
显示结果如下：
D =
0
L =
0.3636 0.5000 1.0000
0.2727 1.0000 0
1.0000 0 0
U =
11.0000 3.0000 0
0 -1.8182 2.0000
0 0 0.0000
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 2.018587e-017.
> In D:\Matlab\pujun\lx0720.m at line 4
X =
1.0e+016 *
-0.4053
1.4862
1.3511
说明：结果中的警告是由于系数行列式为零产生的。可以通过A*X验证其正确性。
（2）Cholesky分解
若A为对称正定矩阵，则Cholesky分解可将矩阵A分解成上三角矩阵和其转置的乘积，即： 其中R为上三角阵。
方程 A*X=b 变成
所以
R = chol(X)
（3）QR分解
对于任何长方矩阵A，都可以进行QR分解，其中Q为正交矩阵，R为上三角矩阵的初等变换形式，即：A=QR
方程 A*X=b 变形成 QRX=b
所以 X=R\(Q\b)
上例中 [Q, R]=qr(A)
X=R\(Q\B)
说明：这三种分解，在求解大型方程组时很有用。其优点是运算速度快、可以节省磁盘空间、节省内存。
5.2 求线性齐次方程组的通解
在Matlab中，函数null用来求解零空间，即满足A·X=0的解空间，实际上是求出解空间的一组基（基础解系）。
格式 z = null % z的列向量为方程组的正交规范基，满足。
% z的列向量是方程AX=0的有理基
例 求解方程组的通解：
解：
>>A=[1 2 2 1;2 1 -2 -2;1 -1 -4 -3];
>>format rat %指定有理式格式输出
>>B=null(A,'r') %求解空间的有理基
运行后显示结果如下：
B =
2 5/3
-2 -4/3
1 0
0 1
或通过行最简行得到基：
>> B=rref(A)
B =
1.0000 0 -2.0000 -1.6667
0 1.0000 2.0000 1.3333
0 0 0 0
即可写出其基础解系（与上面结果一致）。
写出通解：
syms k1 k2
X=k1*B(:,1)+k2*B(:,2) %写出方程组的通解
pretty(X) %让通解表达式更加精美
运行后结果如下：
X =
[ 2*k1+5/3*k2]
[ -2*k1-4/3*k2]
[ k1]
[ k2]
% 下面是其简化形式
[2k1 + 5/3k2 ]
[ ]
[-2k1 - 4/3k2]
[ ]
[ k1 ]
[ ]
[ k2 ]
5.3 求非齐次线性方程组的通解
非齐次线性方程组需要先判断方程组是否有解，若有解，再去求通解。
因此，步骤为：
第一步：判断AX=b是否有解，若有解则进行第二步
第二步：求AX=b的一个特解
第三步：求AX=0的通解
第四步：AX=b的通解= AX=0的通解+AX=b的一个特解。
例 求解方程组
解：见ex06.m
例 求解方程组的通解：（见ex07.m）
解法一：在Matlab编辑器中建立M文件如下：
A=[1 1 -3 -1;3 -1 -3 4;1 5 -9 -8];
b=[1 4 0]';
B=[A b];
n=4;
R_A=rank(A)
R_B=rank(B)
format rat
if R_A==R_B&R_A==n
X=A\b
elseif R_A==R_B&R_Aend
运行后结果显示为：
R_A =
2
R_B =
2
Warning: Rank deficient, rank = 2 tol = 8.8373e-015.
> In ex07 at 13
X =
0
0
-8/15
3/5
C =
3/2 -3/4
3/2 7/4
1 0
0 1
所以原方程组的通解为X=k1+k2+
解法二：用rref求解
A=[1 1 -3 -1;3 -1 -3 4;1 5 -9 -8];
b=[1 4 0]';
B=[A b];
C=rref(B) %求增广矩阵的行最简形，可得最简同解方程组。
运行后结果显示为：
C =
1 0 -3/2 3/4 5/4
0 1 -3/2 -7/4 -1/4
0 0 0 0 0
对应齐次方程组的基础解系为：， 非齐次方程组的特解为：所以，原方程组的通解为：X=k1+k2+。
5.4 非线性方程的解
非线性方程的标准形式为f(x)=0
函数 fzero
格式 x = fzero (fun,x0) %用fun定义表达式f(x)，x0为初始解。
x = fzero (fun,x0,options)
[x,fval] = fzero(…) %fval=f(x)
[x,fval,exitflag] = fzero(…)
[x,fval,exitflag,output] = fzero(…)
说明 该函数采用数值解求方程f(x)=0的根。
例 求的根
解： fun='x^3-2x-5';
z=fzero(fun,2) %初始估计值为2
5.5 非线性方程组的解
非线性方程组的标准形式为：F(x) = 0
其中：x为向量，F(x)为函数向量。
函数 fsolve
格式 x = fsolve(fun,x0) %用fun定义向量函数，其定义方式为：先定义方程函数function F = myfun (x)。
F =[表达式1；表达式2；…表达式m] %保存为myfun.m，并用下面方式调用：x = fsolve(@myfun,x0)，x0为初始估计值。
x = fsolve(fun,x0,options)
[x,fval] = fsolve(…) %fval=F(x)，即函数值向量
[x,fval,exitflag] = fsolve(…)
[x,fval,exitflag,output] = fsolve(…)
例 求下列方程组的根
解：化为标准形式
设初值点为x0=[-5 -5]。
先建立方程函数文件，并保存为myfun.m：
function F = myfun(x)
F = [2*x(1) - x(2) - exp(-x(1));
-x(1) + 2*x(2) - exp(-x(2))];
然后调用优化程序
x0 = [-5; -5]; % 初始点
options=optimset('Display','iter'); % 显示输出信息
[x,fval] = fsolve(@myfun,x0,options)

6. 常微分方程求解
求微分方程的解析解：
r = dsolve('eq1,eq2,...', 'cond1,cond2,...', 'v')
r = dsolve('eq1','eq2',...,'cond1','cond2',...,'v')
dsolve('eq1,eq2,...', 'cond1,cond2,...', 'v')
注意：在表达微分方程时，用字母D表示求微分，D2、D3等表示求高阶微分。任何D后所跟的字母为因变量，自变量可以指定或由系统规则选定为确省。
例如，微分方程 应表达为：D2y=0.
dsolve('Dy=x','x')
例:求微分方程的特解.
y=dsolve('D2y+4*Dy+29*y=0' ,'x')
y=dsolve('D2y+4*Dy+29*y=0','y(0)=0,Dy(0)=15','x')
微分方程数值解法：
[t,x]=ode23(‘xfun’,[t0 tf],y0,tol)
[t,x]=ode45(‘xfun’, [t0 tf],y0,tol)
其中：xfun为以待解方程或方程组写成的M函数
[t0 tf]:自变量初值和终值
yo:函数的初值
tol:设置容许误差（相对误差和绝对误差）
求解过程：
options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-5]);
[T,Y] = ode45(@rigid,[0 12],[0 1 1],options);
plot(T,Y(:,1),'-',T,Y(:,2),'-.',T,Y(:,3),'.')
四、数据处理
插值和拟合都是数据优化的一种方法，当实验数据不够多时经常需要用到这种方法来画图。在matlab中都有特定的函数来完成这些功能。这两种方法的确别在于：
当测量值是准确的，没有误差时，一般用插值；
当测量值与真实值有误差时，一般用数据拟合。
1. 插值
对于一维曲线的插值，一般用到的函数yi=interp1(X,Y,xi,method) ，其中method包括nearst，linear，spline，cubic。
命令1 interp1
功能 一维数据插值。该命令对数据点之间计算内插值。它找出一元函数f(x)在中间点的数值。其中函数f(x)由所给数据决定。各个参量之间的关系如下图。

图1 数据点与插值点关系示意图
格式 yi = interp1(x,Y,xi)
%返回插值向量yi，每一元素对应于参量xi，同时由向量x与Y的内插值决定。
yi = interp1(Y,xi)
%假定x=1:N，其中N为向量Y的长度，或者为矩阵Y的行数。
yi = interp1(x,Y,xi,method) %用指定的算法计算插值：
’nearest’：最近邻点插值，直接完成计算；
’linear’：线性插值（缺省方式），直接完成计算；
’spline’：三次样条函数插值；
’pchip’：分段三次Hermite插值；
’cubic’：与’pchip’操作相同；
’v5cubic’：在MATLAB 5.0中的三次插值。
对于超出x范围的xi的分量，使用方法’nearest’、’linear’、’v5cubic’的插值算法，相应地将返回NaN。对其他的方法，interp1将对超出的分量执行外插值算法。
例 x = 0:10; y = x.*sin(x);
xx = 0:.25:10; yy = interp1(x,y,xx);
plot(x,y,'kd',xx,yy)
例 year = 1900:10:2010;
product = [75.995 91.972 105.711 123.203 131.669 150.697 179.323 203.212 226.505 249.633 256.344 267.893 ];
p1995 = interp1(year,product,1995)
x = 1900:1:2010;
y = interp1(year,product,x,'pchip');
plot(year,product,'o',x,y)

对于二维曲面的插值，一般用到的函数zi=interp2(X,Y,Z,xi,yi,method)，其中method也和上面一样，常用的是cubic。
命令2 interp2
功能 二维数据内插值
格式 ZI = interp2(X,Y,Z,XI,YI)
%返回矩阵ZI，其元素包含对应于参量XI与YI（可以是向量、或同型矩阵）的元素，即Zi(i,j)←[Xi(i,j),yi(i,j)]。
ZI = interp2(Z,XI,YI)
%缺省地，X=1:n、Y=1:m，其中[m,n]=size(Z)。再按第一种情形进行计算。
ZI = interp2(Z,n)
%作n次递归计算，在Z的每两个元素之间插入它们的二维插值，这样，Z的阶数将不断增加。interp2(Z)等价于interp2(z,1)。
ZI = interp2(X,Y,Z,XI,YI,method) %用指定的算法method计算二维插值：
’linear’：双线性插值算法（缺省算法）；
’nearest’：最临近插值