matlab数值计算与符号计算.ppt

第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。 例如，多项式 在MATLAB中用向量p=1 -12 0 25 116表 示。 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 。 (3) 多项式的导函数 对多项式求导数的函数有 k=polyder(p)，返回多项式p的导函数； k=polyder(a, b)，返回多项式a与b的乘积 的导函数； q, d=polyder(b, a)，返回多项式b整除a 的导函数，其分子多项式返回给q，分母 多项式返回给d。 (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) 例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) 返回 A= 16 16 16 16 16 15 -140 -563 16 -140 -2549 -12089 16 -563 -12089 -43779 B= 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 (5) 多项式的根 使用函数roots可以求出多项式等于0的根 ，根用列向量表示，其调用格式为 r=roots(p) 若已知多项式等于0的根，函数poly可以 求出相应多项式，调用格式为 p=poly(r) 例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 3曲线拟合 在MATLAB中用polyfit函数来求得最小二乘拟合 多项式的系数，再用polyval函数按所得的多项 式计算所给出的点上的函数近似值。 polyfit函数的调用格式为 P, S=polyfit(X, Y, m) 函数根据采样点X和采样点函数值Y，产生一个 m次多项式P及其在采样点的误差向量S。其中 X, Y是两个等长的向量，P是一个长度为m + 1 的向量，P的元素为多项式系数。 例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 绘出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”错误。 例4.4 向量t和p表示从19001990年的每隔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年的人口，并利用 插值估计从19002000年每年的人口数。 首先在命令窗口输入插值命令interp1(t, p, 1975) ，估计1975年的人口，返回结果如下： ans= 214.8585 再利用插值估计19002000年每年的人口数， 并利用画图命令得出曲线分布。在MATLAB命 令窗口输入如下命令语句： x=1900:1:2000; y=interp1(t, p, x, spline);%样条插值 plot(t, p, :o, x, y, -r)%带圆圈标记的 虚线(:o)为普查数据，红实线 %(-r)为插值数据 所得曲线如图4.2所示。 图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”错误。 例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.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年所获得的工资。 MATLAB程序如下： w=interp2(service, years, wage, 15, 1975) 程序运行结果如下： w= 190.6287 4.2.1 数值微分 DX=diff(X)，计算向量X的向前差分， 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.2 数值微积分 例4.7 设x由0, 2间均匀分布的10个点组成， 求sinx的13阶差分。 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) 程序运行结果如下： 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.2 数值积分 quad(filename, a, b, tol, trace) quadl(filename, a, b, tol, trace) 其中filename是被积函数名。a和b分别是 定积分的下限和上限。tol用来控制积分精 度，默认时为10-6。trace控制是否展现积 分过程，若取非0则展现积分过程，若取0 则不展现，默认时取trace=0。 例4.9 用两种方法求 在MATLAB命令窗口输入 F=inline(1./(x.3-2*x-5); Q=quad(F, 0, 2) 或 Q=quadl(F, 0, 2) 返回 Q= -0.4605 或者采用函数句柄 Q=quad(myfun, 0, 2) 也可以返回 Q= -0.4605 这里myfun.m为一个M文件： function y=myfun(x) y=1./(x.3-2*x-5); 二重定积分 I=dblquad(f, a, b, c, d, tol, trace) 该函数求f在a, b c, d区域上的二重定 积分。参数tol, trace的用法与函数quad完 全相同。 例4.11 计算二重定积分 首先建立一个函数文件fxy.m， function f=fxy(x, y) global ki; ki=ki+1; %ki用于统计被积函数的调用次 数 f=exp(-x.2/2).*sin(x.2+y); 再调用dblquad函数求解。 global ki;ki=0; I=dblquad(fxy, -2, 2, -1, 1) ki 程序运行结果如下： I= 1.57449318974494 ki= 1038 4.3 线性方程组求解 对于线性方程组Ax=b，可以利用左除运 算符“”求解： x=Ab 即x=inv(A)*b 如果矩阵A是奇异的或接近奇异的，则 MATLAB会给出警告信息。 MATLAB提供的lu函数用于对矩阵进行LU分解 ，其调用格式如下。 L, U=lu(X)，产生一个上三角阵U和一个变换形 式的下三角阵L(行交换)，使之满足X=LU。注意 ，这里的矩阵X必须是方阵。 L, U, P=lu(X)，产生一个上三角阵U和一个下 三角阵L以及一个置换矩阵P，使之满足PX=LU 。当然矩阵X同样必须是方阵。 实现LU分解后，线性方程组Ax=b的解可表示为 x=U(Lb)或x=U(LPb)，这样可以大大提高运 算速度。 MATLAB的函数qr可用于对矩阵进行QR 分解，其调用格式如下。 Q, R=qr(X)，产生一个正交矩阵Q和一个 上三角矩阵R，使之满足X=QR。 Q, R, E=qr(X)，产生一个正交矩阵Q、 一个上三角矩阵R以及一个置换矩阵E， 使之满足XE=QR。 实现QR分解后，线性方程组Ax=b的解可 表示为x=R(Qb)或x=E(R(Qb)。 MATLAB函数chol(X)用于对矩阵X进行 Cholesky分解，其调用格式如下。 R=chol(X)，产生一个上三角阵R，使RR=X。 若X为非对称正定的，则输出一个出错信息。 R, p=chol(X)，这个命令格式将不输出出错信 息。当X为对称正定的，则p=0，R与上述格式 得到的结果相同；否则p为一个正整数。如果X 为满秩矩阵，则R为一个阶数为q=p-1的上三角 阵，且满足RR=X(1:q, 1:q)。 实现Cholesky分解后，线性方程组Ax=b变成 RRx=b，所以x=R(Rb)。 例4.12 用直接法求解下列线性方程 组。 MATLAB程序如下： A=1/2 1/3 1/4;1/3 1/4 1/5;1/4 1/5 1/6 B=0.95 0.67 0.52 X=AB或X=inv(A)*B 程序运行结果如下： X= 1.2000 0.6000 0.6000 例4.13 分别用LU分解、QR分解和Cholesky 分解求解例4.12中的线性方程组。 MATLAB程序如下： A=1/2 1/3 1/4;1/3 1/4 1/5;1/4 1/5 1/6 B=0.95 0.67 0.52 L, U=lu(A); X=U(LB) %LU分解求解 Q, R=qr(A); X=R(QB) %QR分解求解 R=chol(A) X=R(RB) %Cholesky分解求解 所得结果均为 X= 1.2000 0.6000 0.6000 Jacobi迭代法的MATLAB函数文件Jacobi.m如下： function y, n=jacobi(A, b, x0, eps) if nargin=3 eps=1.0e-6; elseif nargin=eps x0=y; y=B*x0+f; n=n+1; end 例4.14 用Jacobi迭代法求解线性方程组。设 迭代初值为0，迭代精度为10-6。 在命令中调用函数文件Jacobi.m，命令如 下： A=9, -1, 1;-1, 10, -2; -2, 1, 10; b=10, 7, 6; x, n=jacobi(A, b, 0, 0, 0, 1.0e-6) MATLAB返回 x= 1.1365 0.9599 0.7313 n= 10 输入x=inv(A)*b，所得结果相同。 4.4 常微分方程的数值求解 MATLAB提供了求常微分方程数值解的函数， 一般调用格式如下： t, y=ode23(fname, tspan, y0), 二阶、三阶龙 格-库塔法 t, y=ode45(fname, tspan, y0)，四阶、五阶龙 格-库塔法 其中，fname是定义f(t, y)的函数文件名，该函 数文件必须返回一个列向量；tspan形式为t0, tf；表示求解区间；y0是初始状态列向量；t和y 分别给出时间向量和相应的状态向量。 例4.17 考虑著名的Rossler微分方程组 选定a = b = 0.2，c = 5.7，且x(0) = y(0) = z(0)=0, 求解该微分方程。 引入新状态变量 其矩阵形式为 编写M函数如下： function dx=rossler(t, x, flag, a, b, c) %虽然不显含时间，但还应该写出占位 dx=-x(2)-x(3);x(1)+a*x(2);b+(x(1)- c)*x(3); %对应方程 保存为rossler.m文件，在MATLB命令窗口输入如 下语句： a=0.2;b=0.2;c=5.7; %从函数外部定义参数变量 x0=0 0 0; %微分方程的初值 t, y=ode45(rossler, 0, 100, x0, , a, b, c) %求解微分方程 plot(t, y) %绘制各个状态变量的时间响应 figure plot3(y(:, 1), y(:, 2), y(:, 3) 上面的命令直接得出该微分方程在内的数 值解，并绘制出各个状态变量和时间之间 的关系曲线和相空间曲线，如图4.5(a)和 图4.5(b)所示。 图4.5 Rossler方程的数值解表示 4.5 MATLAB符号计算 (1) sym函数 sym函数用来建立单个符号量，其调用格 式为 符号变量名=sym(符号字符串) 符号字符串可以是常量、变量、函数或表 达式。 (2) syms函数 syms函数的一般调用格式为 syms var1 var2 varn 函数定义符号变量var1, var2, , varn等。用这 种格式定义符号变量时不要在变量名上加字符 分界符()，变量间用空格而不要用逗号分隔。 syms x beta等同于 x=sym(x); beta=sym(beta); 建立符号表达式有以下两种方法。 (1) 用sym函数建立符号表达式 命令语句y1=sym(1/sqrt(2*x) 返回y1=1/sqrt(2*x) 命令语句M=sym(a, b;c, d) 返回M= a, b c, d (2) 使用已经定义的符号变量 组成符号表达式 命令语句syms x y v=3*x2-5*y+2*x*y+6 返回v= 3*x2-5*y+2*x*y+6 符号表达式的加、减、乘、除运算可分别 由函数symadd，symsub，symmul和 symdiv来实现，幂运算可以由sympow来 实现。此外，和数值运算一样，也可以用 +、-、*、/、运算实现符号运算。 n, d=numden(s) 该函数提取符号表达式s的分子和分母， 分别将它们存放在n与d中。例如： 命令语句n, d=numden(sym(4/5) 返回 n=4 and d=5 factor(S)，对S分解因式，S是符号表达式 或符号矩阵。 expand(S)，对S进行展开，S是符号表达 式或符号矩阵。 collect(S)，对S合并同类项，S是符号表 达式或符号矩阵。 collect(S, v)，对S按变量v合并同类项，S 是符号表达式或符号矩阵。 例如： 命令语句f=factor(123) 返回f= 3 41 命令语句expand(x-2)*(x-4) 返回x2-6*x+8 命令语句expand(cos(x+y) 返回cos(x)*cos(y)-sin(x)*sin(y) 命令语句expand(exp(a+b)2) 返回exp(a2)*exp(a*b)2*exp(b2) simplify(S)，应用函数规则对S进行化简。 simple(S)，调用MATLAB的其他函数对表达式 进行综合化简，并显示化简过程。 例如：命令语句simplify(sin(x)2 + cos(x)2) 返回1 命令语句simplify(exp(c*log(sqrt(a+b) 返回(a+b)(1/2c) 命令语句S=(x2+5*x+6)/(x+2), sqrt(16); R=simplify(S) 返回R=x+3, 4 利用函数sym可以将数值表达式变换成它 的符号表达式。例如： 命令语句sym(1.5) 返回 ans= 3/2 命令语句sym(3.14) 返回ans= 157/50 函数numeric或eval可以将符号表达式变换成数值表达 式。N=numeric(S)等同于N=double(sym(S)，例如： 命令语句A=magic(4); A(:, :, 2)=A; d1, d2, d3=eval(size(A) 返回d1= 4 d2= 4 d3= 2 transpose(S)，返回S矩阵的转置矩阵。 determ(S)，返回S矩阵的行列式值。 colspace(S)，返回S矩阵列空间的基。 ezplot(f)，绘制函数f(x)在默认区间-2x2 内的图形； ezplot(f, min, max)，绘制函数f(x)在指定区间 min, max内的图形。该函数打开标签为Figure No.1的图形窗口，并显示图形。如果已经存在 图形窗口，该函数在标签数最大的窗口中显示 图形； ezplot(f, xmin xmax, fign)，在指定的窗口fign 中绘制函数f(x)在指定区间min, max内的图形 。 对于隐函数，ezplot函数的调用格式有： ezplot(f)，绘制函数f(x)在区间-2x2 内的图形； ezplot(f, xmin, xmax, ymin, ymax)，绘 制函数f(x, y)在指定区间xminxxmax ，yminyymax的图形； ezplot(f, min, max)，绘制函数f(x, y)在指 定区间minxmax，minymax的图 形。 对于参数方程，ezplot函数的调用格式有 ： ezplot(x, y)，绘制参数方程x=x(t)，y=y(t) 在默认区间0t2内的图形； ezplot(x, y, tmin, tmax)，绘制参数方程 x=x(t)，y=y(t)在指定区间tinttmax的 图形。 ezplot3函数用于绘制三维参数曲线。该函 数的调用格式有： ezplot3(x, y, z)，在默认区间0t2内 绘制参数方程x=x(t)，y=y(t)，z=z(t)的图 形； ezplot3(x, y, z, tmin, tmax)，在区间tmin ttmax内绘制参数方程x=x(t)，y=y(t) ，z=z(t)的图形； ezplot3(., animate)，生成空间曲线的动 态轨迹。 ezmesh，ezsurf函数分别用于绘制三维网格曲 面图和三维表面图。这两个函数的用法相同， ezmesh函数的调用格式有： ezmesh(f)，绘制函数f(x, y)在 默认区间-2x 2，-2y2内的图形； ezmesh(f, domain)，在指定区域domain绘制函 数f(x, y)的图形； ezmesh(x, y, z)，绘制参数方程x=x(s，t)， y=y(s，t)，z=y(s，t)在默认区间-2s2，- 2t2内的图形 ezmesh(x, y, z, smin, smax, tmin, tmax)或 ezmesh(x, y, z, min, max)，在指定区域绘制 三维参数方程的图形。 ezmeshc，ezsurfc两个函数用于在绘制三维曲 面的同时绘制等值线。 ezmeshc(f)，绘制二元函数f(x,y)在默认区间-2 x2，-2y2的图形。 ezmeshc(f, domain)，绘制函数f(x,y)在指定区 域的图形，绘图区域由domain指定，其中 domain为4 1数组或者2 1数组，如xmin, xmax, ymin, ymax表示minxmax，miny max，min, max表示minxmax, miny max。 ezmeshc(x, y, z)，绘制参数方程x=x(s, t)， y=y(s, t)和z=z(s, t)在默认区间-2s2，-2 t2的图形。 ezmeshc(x, y, z, smin, smax, tmin, tmax)， ezmeshc(x, y, z, min, max)，绘制参数方程在 指定区域的图形。 ezmeshc(., n)，指定绘图的网格数，默认值为 60。 ezmeshc(., circ)，在以指定区域中心为中心 的圆盘上绘制图形 在MATLAB中，用于绘制符号函数等值线的函数有 ezcontour和ezcontourf，这两个函数分别用于绘制等值 线和带有区域填充的等值线。其调用格式有： ezcontour(f)，绘制符号二元函数f(x, y)在默认区间-2 x2，-2y2内的等值线图； ezcontour(f, domain)，绘制符号二元函数f(x, y)在指定 区域的等值线图； ezcontour(., n)，绘制等值线图，并指定等值线的数目 。 例如命令ezplot(x2-y4)绘制的图形曲线如图4.7所示 。 例如命令ezmeshc(y/(1+x2+y2), -5, 5, -2*pi, 2*pi)绘 制的图形如图4.8所示。 limit(f, x, a)，计算当变量x趋近于常数a时，f(x) 函数的极限值。 limit(f, a)，求符号函数f(x)的极限值，符号函数 f(x)的变量为函数findsym(f)确定的默认自变量， 即变量x趋近于a。 limit(f)，系统默认变量趋近于0，即a=0的极限 。 limit(f, x, a, right)，变量x从右边趋近于a时符号 函数f(x)的极限值。 limit(f, x, a, left)，变量x从左边趋近于a时符号 函数的极限值。 diff(f, x, n) diff函数求函数f对变量x的n阶导数。参数x 的用法同求极限函数limit，可以默认，默 认值与limit相同，n的默认值是1。 int(f, x)，求函数f对变量x的不定积分。参 数x可以默认，默认原则与diff函数相同。 int(f, v)，以v为自变量，对被积函数符号 表达式f求不定积分。 int(f, v, a, b)，求被积函数f在区间a, b上 的定积分。a和b可以是两个具体的数，也 可以是一个符号表达式。 例4.20 求椭球