数值计算与符号计算

• • • • •
例4.9 用两种方法求 I 0 在MATLAB命令窗口输入 F=inline('1./(x.^3-2*x-5)'); Q=quad(F, 0, 2) 或 Q=quadl(F, 0, 2)
2
1 dx 3 x 2x 5
• • • • • • • •
返回 Q= -0.4605 或者采用函数句柄 Q=quad(@myfun, 0, 2) 也可以返回 Q= -0.4605
• • • • •
MATLAB程序如下： w=interp2(service, years, wage, 15, 1975) 程序运行结果如下： w= 190.6287
4.2 数值微积分
• 4.2.1 数值微分 • DX=diff(X)，计算向ຫໍສະໝຸດ BaiduX的向前差分， DX(i)=X(i+1)-X(i)，i=1, 2, …, n-1。 • DX=diff(X, n)，计算X的n阶向前差分，例 如diff(X, 2)=diff(diff(X))。 • DX=diff(A, n, dim)，计算矩阵A的n阶差分， dim=1时(默认状态)，按列计算差分； dim=2时，按行计算差分。
• 例4.4 向量t和p表示从1900～1990年的每隔10年 的美国人口普查数据： • t=1900:10:1990; • p=[75.995 91.972 105.711 123.203 131.669... • 150.697 179.323 203.212 226.505 249.633]; • 根据人口普查数据估计1975年的人口，并利用插 值估计从1900～2000年每年的人口数。
• 例4.5 利用插值运算对峰值函数peaks插入更多的 栅格。 • MATLAB程序如下： [X, Y]=meshgrid(-3:.25:3); • Z=peaks(X, Y); • [XI, YI]=meshgrid(-3:.125:3); • ZI=interp2(X, Y, Z, XI, YI); • mesh(X, Y, Z), hold, • mesh(XI, YI, ZI+15) • hold off • axis([-3 3 -3 3 -5 20]) • 程序运行结果如图4.3所示。
• 例4.7 设x由[0, 2]间均匀分布的10个点组成，求 sinx的1～3阶差分。 • MATLAB程序如下： • 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)
• • • • • • • • • • • • • • •
例4.2 求多项式x4 + 8x3 - 10的根。 命令如下： A=[1, 8, 0, 0, -10]; x=roots(A) 返回 x= -8.0194 1.0344 -0.5075 + 0.9736i -0.5075 - 0.9736i 再输入 p=poly(x) 返回 p= 1.0000 8.0000 -0.0000 -0.0000 -10.0000
例4.3 用一个6次多项式在区间[0, 2]内 逼近函数
• • • • MATLAB程序如下： x=linspace(0, 2*pi, 50); y=sin(x); P=polyfit(x, y, 6) %得到6次多项式的系 数和误差 • 程序运行结果如下： • P= • 0.0000 -0.0056 0.0874 -0.3946 0.2685 0.8797 0.0102
• 返回 • A= • 16 • 16 • 16 • 16 • B= • 0 • 0 • 0 • 0
16 15 -140 -563
0 0 0 0
16 -140 -2549 -12089
0 0 0 0
16 -563 -12089 -43779
0 0 0 0
• (5) 多项式的根 • 使用函数roots可以求出多项式等于0的根， 根用列向量表示，其调用格式为 • r=roots(p) • 若已知多项式等于0的根，函数poly可以求 出相应多项式，调用格式为 • p=poly(r)
3．曲线拟合
• 在MATLAB中用polyfit函数来求得最小二乘拟合多 项式的系数，再用polyval函数按所得的多项式计 算所给出的点上的函数近似值。 • polyfit函数的调用格式为 • [P, S]=polyfit(X, Y, m) • 函数根据采样点X和采样点函数值Y，产生一个m 次多项式P及其在采样点的误差向量S。其中X, Y 是两个等长的向量，P是一个长度为m + 1的向量， P的元素为多项式系数。
4.2.2 数值积分
• quad(filename, a, b, tol, trace) • quadl(filename, a, b, tol, trace) • 其中filename是被积函数名。a和b分别是定 积分的下限和上限。tol用来控制积分精度， 默认时为10-6。trace控制是否展现积分过 程，若取非0则展现积分过程，若取0则不 展现，默认时取trace=0。
• • • •
(3) 多项式的导函数 对多项式求导数的函数有 k=polyder(p)，返回多项式p的导函数； k=polyder(a, b)，返回多项式a与b的乘积的 导函数； • [q, d]=polyder(b, a)，返回多项式b整除a的 导函数，其分子多项式返回给q，分母多项 式返回给d。
• 程序运行结果如下： • X= • 0 0.6981 1.3963 2.0944 2.7925 3.4907 4.1888 4.8869 5.5851 6.2832 • Y= • 0 0.6428 0.9848 0.8660 0.3420 -0.3420 -0.8660 -0.9848 0.6428 -0.0000 • DY= • 0.6428 0.3420 -0.1188 -0.5240 -0.6840 -0.5240 -0.1188 0.3420 0.6428 • D2Y= • -0.3008 -0.4608 -0.4052 -0.1600 0.1600 0.4052 0.4608 0.3008 • D3Y= • -0.1600 0.0556 0.2452 0.3201 0.2452 0.0556 -0.1600
图4.2 用一维数据插值得到的美国100年人口分布图
• (2) 二维数据插值 • 在MATLAB中，提供了解决二维插值问题的函数 interp2，其调用格式为 • Z1=interp2(X, Y, Z, X1, Y1, 'method') • 其中，X, Y是两个向量，分别描述两个参数的采 样点；Z是与参数采样点对应的函数值；X1, Y1是 两个向量或标量，描述欲插值的点；Z1是根据相 应的插值方法得到的插值结果。method的取值与 一维插值函数相同。X, Y, Z也可以是矩阵形式。 • 同样，X1, Y1的取值范围不能超出X, Y的给定范 围，否则，会给出“NaN”错误。
• 首先在命令窗口输入插值命令interp1(t, p, 1975)， 估计1975年的人口，返回结果如下： • ans= • 214.8585 • 再利用插值估计1900～2000年每年的人口数，并 利用画图命令得出曲线分布。在MATLAB命令窗 口输入如下命令语句： • x=1900:1:2000; • y=interp1(t, p, x, 'spline'); %样条插值 • plot(t, p, ':o', x, y, '-r') %带圆圈标记的虚线 (:o)为普查数据，红实线 • %(-r)为插值数据 • 所得曲线如图4.2所示。
• 例4.1 多项式P = x4-29x3 + 72x2 - 29x + 1， 以4阶pascal矩阵为自变量分别用polyval和 polyvalm计算该多项式的值。 • 在命令窗口输入如下命令： • p= [1 -29 72 -29 1]; • X=pascal(4); • A=polyval(p, X), B=polyvalm(p, X)
• 图4.3 用二维数据插值得到的peaks函数曲线图
• • • • • • • • •
例4.6 给定雇员数据如下： years=1950:10:1990; %工作年份 service=10:10:30; %服役时间，即在职时间 wage=[150.697 199.592 187.625 179.323 195.072 250.287 203.212 179.092 322.767 226.505 153.706 426.730 249.633 120.281 598.243]; %工资 利用二维数据插值计算一个雇员在工作15年后在 1975年所获得的工资。
• (4) 多项式求值 • MATLAB中提供了两种求多项式值的函数。 • y=polyval(p, x)，代数多项式函数求值，若x为一 数值，则求多项式在该点的值；若x为向量或矩阵， 则对向量或矩阵中的每个元素求其多项式的值。 • Y=polyvalm(p, x)，矩阵多项式求值，要求x为方 阵。设A为方阵，p代表多项式x3-5x2+8，那么 polyvalm(p, A)的含义是 • A*A*A-5*A*A+8*eye(size(A)) • 而polyval(p, A)的含义是 • A.*A.*A-5*A.*A+8*ones(size(A))
2．多项式的运算
• (1) 多项式的加减运算 • 多项式的加减运算就是其所对应的系数向 量的加减运算。相加减的多项式必须表示 成相同的次数，如果次数不同，应该把低 次的多项式不足的高次项用0补足。
• (2) 多项式的乘除运算 • 命令w=conv(u, v)表示多项式u和v相乘，例如 在MATLAB中输入 • u=[1 2 3 4]，v=[10 20 30]，c=conv(u, v) • 返回 • c= • 10 40 100 160 170 120 • conv指令可以嵌套使用，如conv(conv(a, b), c)。 • 命令[q, r]=deconv(v, u)表示u整除v。 向量q表 示商，向量r表示余， 即有v=conv(u, q)+r。
第4章 MATLAB数值计算与符号计算
4.1曲线拟合与插值运算 4.2 数值微积分 4.3 线性方程组求解 4.4 常微分方程的数值求解 4.5 MATLAB符号计算 4.6 级数 4.7 实验五 数值工具箱与符号工具箱的应用
4.1 曲线拟合与插值运算
• 1．多项式的建立与表示方法 • 在MATLAB中，n次多项式用一个长度为n + 1的行向量表示，其元素为多项式的系数，按 降幂排列，缺少的幂次项系数为0。 • 例如，多项式 x 4 12 x3 0 x 2 25 x 116 • 在MATLAB中用向量p=[1 -12 0 25 116]表示。
• 绘出sinx和多项式P(x)在给定区间的函数曲 线，如图4.1所示。
图4.1 用6次多项式对正弦函数进行拟合
4．数据插值
• (1) 一维数据插值 • 在MATLAB中，实现一维数据插值的函数是interp1，被 插值函数是一个单变量函数，其调用格式为 • yi=interp1(x, Y, xi, method) • 函数根据x, y的值，计算函数在xi处的值。x, y是两个等 长的已知向量，分别描述采样点和样本值，xi是一个向 量或标量，描述欲插值的点，yi是一个与xi等长的插值 结果。method是插值方法，允许的取值有'linear'、 'nearest'、'cubic'、'spline'，分别表示线性插值、最近 点插值、3次多项式插值、3次样条插值。 • 注意：xi的取值范围不能超出x的给定范围，否则，会给 出“NaN”错误。