第3章 MATLAB数值计算 2017/9/9

第3章 MATLAB数值计算 3.1 多项式 3.2 插值和拟合 3.3 数值微积分 3.4 线性方程组的数值解 3.5 稀疏矩阵 3.1 多项式 3.2 插值和拟合 3.3 数值微积分 3.4 线性方程组的数值解 3.5 稀疏矩阵 3.6 常微分方程的数值解

3.1 多项式 3.1.1 多项式的表达和创建 表示成向量的形式，系数按降序排列 例如 x = [1 3 -15 -2 9] y = [1 0 0 0 1]

3.1.2 多项式的四则运算 多项式相加减没有专门的函数，可以自己定义。 多项式相乘函数 conv 的语法为 c=conv(a,b)，其 【例 3.4】 完成两个同阶次多项式： 的相乘运算。 >> a=[1 2 3 4]; >> b=[1 4 9 16]; >> e = conv(a,b) e = 1 6 20 50 75 84 64

多项式相除函数 deconv 的语法为 [q,r]=deconv(a,b)，其中 【例 3.6】 利用例 3.4 中的数据。 >> [f, r] = deconv(e,b) f = 1 2 3 4 r = 0 0 0 0 0 0 0

3.1.3 多项式求值和求根运算 1. 多项式求值 语法格式为 y = polyval(p,x) 1. 多项式求值 语法格式为 y = polyval(p,x) 其中 p 代表多项式各阶系数向量，x 为要求值的点。当 x 表示矩阵时，需用 y=polyvalm(p,x)来计算相应的值。 【例 3.8】 利用 polyval 函数找出 在 s=3 处的值： >>p=[1 2 -12 -1 7]; >>z=polyval(p,3) z = 31

【例 3.9】 利用 polyval找出多项式 在[-1,4]间均匀分布的 5个离散点的值。 >> x=linspace(-1,4,5) % 在[-1,4]区间产生5个离散点 >> p=[1 4 7 -8]; >> v=polyval(p,x) x = -1.0000 0.2500 1.5000 2.7500 4.0000 v = -12.0000 -5.9844 14.8750 62.2969 148.0000 v 即为多项式在各个离散点上对应的函数值。

【例 3.10】 估计矩阵多项式 P (X) = X^3 – 2X – I在已知矩阵 X 处的值，其中 X=[1 2 1; -1 0 2; 4 1 2]。 >>Y = polyvalm(P,X) Y = 0 -1 5 9 -1 -1 3 8 5 2. 多项式求根 语法为：x=roots(P)，其中 P 为多项式的系数向量，x 也为向量，即 x(1),x(2),…,x(n)分别代表多项式的 n 个根。MATLAB规定：多项式是行向量，根是列向量。

【例 3.11】 求解多项式 的根。 >>roots([1 3 -12 -2 8]) ans = -5.18325528043789 2.17062070347062 -0.83694739215044 0.84958196911772 注意：在上面的程序中，数字格式都设为长(long)型，若改为短(short)型，结果会有差别， 根据需要可执行 MATLAB 窗口的 Fle | Preferences命令进行修改。

3.1.4 多项式的构造 函数 poly2sym来构造多项式 函数 poly来求根对应的多项式的各阶系数 >> T=[1 3 -15 -2 9]; >> poly2sym(T); ans = x^4+3*x^3-15*x^2-2*x+9

3.2 插值和拟合 【例 3.15】 用多项式的根构造上例多项式 。 【例 3.15】 用多项式的根构造上例多项式 。 >> T=[1 3 -15 -2 9]; %多项式的系数向量 >> r=roots(T); %求得多项式的根 >> poly(r) %利用根构造出多项式 ans = 1.0000 3.0000 -15.0000 -2.0000 9.0000 3.2 插值和拟合 3.2.1 多项式插值和拟合 插值 已知 节点 构造函数 使得

拟合 拟合就是要找出一个曲线方程式（多项式拟合就是设 法找一个多项式），使得它与观测数据最为接近，这时 不要求拟合多项式通过全部已知的观测节点。 1．多项式插值函数(interp1) yi = interp1(x,y,xi,method) 对应于插值函数 ，其中 x 和 y是原已知数据的 x、y 值，xi 是要内插的数据点，method是插值方法。

【例 3.16】 取余弦曲线上 11 个点的自变量和函数值点作为已知数据，再选取 41 个自变量点，分别用分段线性插值、三次方程式插值和样条插值3 种方法计算确定插值函数的值。 >> x=0:10; y=cos(x); >> xi=0:.25:10; >> y0=cos(xi); %精确值 >> y1=interp1(x,y,xi); %线性插值结果 >> y2=interp1(x,y,xi,'cubic'); %三次方程式插值结果 >> y3=interp1(x,y,xi,'spline'); %样条插值结果 >> plot(xi,y0,'o',xi,y1,xi,y2,'-.',xi,y3)

2．多项式拟合函数 polyfit p=polyfit(x,y,n) [p,s]=polyfit(x,y,n) 其中x,y为已知的数据组，n 为要拟合的多项式的阶次，向量p 为返回的要拟合的多项式的系数，向量s 为调用函数polyval 获得的错误预估计值。 一般来说，多项式拟合中阶数n 越大，拟合的精度就越高。 函数 polyfit 常和函数 polyval(见 3.1.3 节)结合起来使用，由 polyfit 计算出多项式的各个系数 后，再利用polyval 对输入向量决定的多项式求值。

【例 3.18】 对向量X=[-2.8 -1 0.2 2.1 5.2 6.8]和Y=[3.1 4.6 2.3 1.2 2.3 -1.1]分别进行阶数为3、4、5 的多项式拟合，并画出图形进行比较。 >> p3=polyfit(x, y, 3); % 用不同阶数的多项式拟合x和y >> p4=polyfit(x, y, 4); >> p5=polyfit(x, y, 5); >> xcurve= -3.5:0.1:7.2; % 生成x值 >> p3curve=polyval(p3, xcurve); % 计算在这些x点的多项式值 >> p4curve=polyval(p4, xcurve); >> p5curve=polyval(p5, xcurve); >> plot(xcurve,p3curve,'--',xcurve,p4curve,'-.',xcurve,p5curve,'-',x,y,'*');

3.3 数值微积分 3.3.1 微分和差分 函数diff:计算两个相邻点的差值. 语法 diff(x) —返回x 对预设独立变量的一次微分； 3.3.1 微分和差分 函数diff:计算两个相邻点的差值. 语法 diff(x) —返回x 对预设独立变量的一次微分； diff(x,'t') —返回x 对独立变量t 的一次微分值； diff(x,n) —返回x 对预设独立变量的n 次微分值； diff(x,'t',n) —返回x 对独立变量t 的n 次微分值。 其中x代表一组离散点 。 计算 的数值微分为dy=diff(y)./diff(x)。

【例 3.22】 计算多项式 在 [-4, 5] 区间的微分。 >> x=linspace(-4,5); %产生100个x的离散点 >> p= [1 -3 -8 7 3 -5]; >> f=polyval(p,x); %多项式在100个离散x点上对应的值 >> subplot(2,1,1);plot(x,f) %将多项式函数绘图 >> title('多项式方程') ; >> dfb=diff(f)./diff(x); %注意要分别计算diff(f)和diff(x) >> xd=x(2:length(x)); %注意只有99个df值，而且是对应x2,x3, …, %x100的点 >> subplot(2,1,2);plot(xd,dfb ); %绘多项式的微分图 >> title('多项式方程的微分图') ;

【例 3.23】 对 3 个方程式 利用diff 的4 种语法格式计算微分的示例。 >> S1 = '6*x^3-4*x^2+b*x-5'; %符号表达式(见第5章) >> S2 = 'sin(a)'; >> S3 = '(1 - t^3)/(1 + t^4)'; >> diff(S1) %对预设独立变量x的一次微分值 ans= 18*x^2-8*x+b >> diff(S1,2) %对预设独立变量x的二次微分值 36*x-8 >> diff(S1,'b') %对独立变量b的一次微分值

ans= x >> diff(S2) %对预设独立变量a的一次微分值 cos(a) >> diff(S3) %对预设独立变量t的一次微分值 -3*t^2/(1+t^4)-4*(1-t^3)/(1+t^4)^2*t^3 3.3.2 牛顿-科茨系列数值积分公式 1. 矩形法数值积分 函数cumsum cumsum(x)*h，其中h 为子区间步长，cumsum(x)对应 。

【例 3.25】 利用矩形法计算积分 (该积分的精确值为2)。 >> x=linspace(0,pi,100); %在[0，π]之间取100个离散点 >> y=sin(x); >> T=cumsum(y)*pi/(100-1); % pi/(100-1)表示两个离散点之间的距离 >> I=T(100) %函数在[0，π]之间的矩形积分 I= 1.9998 2. 梯形法数值积分 函数trapz z=trapz(x,y) 表示通过梯形积分法计算 y 对 x 的数值积分。

【例 3.26】 利用梯形法计算积分 。 >> x=linspace(0,pi,100); >> y=sin(x); >> t=trapz(x,y) t = 1.9998 3. 辛普森数值积分 q=quad('f',a,b)

【例 3.28】 用辛普森积分公式求 的积分。 解：方法1 >> quad('1./(x.^3-2*x-5)',0,2) ans = -0.4605 方法2 >> F='1./(x.^3-2*x-5)'; >> quad(F,0,2) ans= 4. 科茨数值积分 函数quadl q = quadl(fun,a,b)

3.4 线性方程组的数值解 【例 3.30】 用科茨积分公式求 的积分。 【例 3.30】 用科茨积分公式求 的积分。 >> quadl('1./(x.^3-2*x-5)',0,2) ans = -0.4605 3.4 线性方程组的数值解 3.4.1 直接法 1. 矩阵相除法 对线性方程组 AX=B X=A\B 2. 消去法 方程的个数和未知数个数不相等，用消去法。将增广矩阵(由[A B]构成)化为简化阶梯形，若系数矩阵的秩不等于增广矩阵的秩，则方程组无解；若两者的秩相等，则方程组有解，方程组的解就是行简化阶梯形所对应的方程组的解。

【例 3.34】 解方程组： >> a=[1 -1 1 -1 1;1 -1 -1 1 0;1 -1 -2 2 -0.5]; %为增广矩阵，由[A B]构成 >> rref(a) ans = 1.0000 -1.0000 0 0 0.5000 0 0 1.0000 -1.0000 0.5000 0 0 0 0 0 由结果看出，x2、x4 为自由未知量，方程组的通解为：x1= x2+0.5，x3= x4+0.5。

3.6 常微分方程的数值解 利用函数dsolve 可求解微分方程(组)的解析解 龙格-库塔方法 3.6 常微分方程的数值解 利用函数dsolve 可求解微分方程(组)的解析解 龙格-库塔方法 最常用的函数是 ode23(二三阶龙格-库塔函数)和 ode45(四五阶龙格-库塔函数) 【例 3.47】 分别用二三阶龙格-库塔法和四五阶龙格-库塔法解常微分方程的初值问题： 取步长 。方程的解析解可用 dsolve求得为 ，将龙格-库塔法与精确值进行比较。

解：先将微分方程写成自定义函数 exam2fun.m function f=exam2fun (x,y) f=-y-x*y.^2; f=f(:); 然后在命令窗口输入以下语句： >> [x1,y1]=ode23('exam2fun',[0:0.1:1],1) >> [x2,y2]=ode45('exam2fun',[0:0.1:1],1)

x1 = 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000 y1 = 1.0000 0.9006 0.8046 0.7144 0.6314 0.5563 0.4892 0.4296 0.3772 0.3312 0.2910 x2 = 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000 y2 = 1.0000 0.9006 0.8046 0.7144 0.6315 0.5563 0.4892 0.4296 0.3772 0.3312 0.2910