matlab第八讲数值计算.ppt

Matlab第八讲 数值计算,赵瑞,数值计算,线性方程组的解 函数的零点 函数的极值点 数值积分 曲线的多项式拟合,一、线性方程组的解,线性方程组AX=b,A为可逆矩阵，求X的命令： X=inv(A)*b %b输入成列向量，X返回列向量 不推荐使用逆矩阵法 2.推荐使用“除法”，精度比逆矩阵法高 X=Ab %A左除b,相当于A-1b,b输入成列向量，X返回列向量 3. X=b/A %b右除A，相当于bA-1, b输入成行向量，X返回行向量 例：A=1,1/2,1/3;1/2,1/3,1/4;1/3,1/4,1/5;b=1;2;3; X1=inv(A)*b X2=Ab b=b; X3=b/A,二、函数的零点,对于任意函数f(x)=0来说，它可能有零点，也可能没有 零点；可能只有一个零点，也可能有多个甚至无数个零 点。因此，很难说出一个通用解法。一般说来，零点的数 值计算过程是：先猜测一个初始零点或该零点所在的区 间；然后通过一些计算，使猜测值精确化。,利用fzero指令求一元函数零点 x,y=fzero(fun,x0) fzero是根据函数是否穿越横轴来决定零点的，因此本程序无法确定函数曲线仅触及横轴和不穿越横轴的零点。例如：y=x2的零点，不穿越横轴。 fun表示函数，数据类型是字符串 x0是表示零点初始猜测。x0可以是标量或二元向量。当x0取标量时，该指令将在它两侧寻找一个与之最靠近的零点；当x0取向量a,b时，该指令将在区间a,b内寻找一个零点。 返回值x表示零点，y表示在这点的函数值 注意：这个命令只能求x0附近的一个零点,例：求函数f(x)=(sin2x)e-0.1x-0.5|x|的零点 (1)首先用字符串表示函数，matlab命令： f=sin(x)2*exp(-0.1*x)-0.5*abs(x) ; （2）作图法观察函数零点分布，matlab命令： x=-10:0.01:10; y=sin(x).2.*exp(-0.1*x)-0.5*abs(x) ; plot(x,y); hold on %接着画,下面画出x轴 xx=-10:0.01:10; yy=zeros(size(xx); plot(xx,yy, r),观察f(x)与x轴的交点，大致可以看出有四个穿越x轴的交点，具体位置不是很清楚，所以希望图像能够局部放大，以便看得更加清楚,（3）利用zoom和ginput指令获得零点的初始近似值，matlab命令： zoom on %在MATLAB指令窗中运行，获局部放大图， %或者直接在图形窗口单击放大镜图标即可,局部放大之后的图形,xx,yy=ginput(5) %在MATLAB指令窗中运行，用鼠标获5个零点猜测值。,zoom off hold off （4）求靠近xx各点的精确零点，matlab命令 x1,y1=fzero(f,xx(1) %求在xx(1)附近的精确零点 x2,y2=fzero(f,xx(2) %求在xx(2)附近的精确零点 x3,y3=fzero(f,xx(3) %求在xx(3)附近的精确零点 x4,y4=fzero(f,xx(4) %求在xx(4)附近的精确零点 x5,y5=fzero(f,xx(5) %求在xx(5)附近的精确零点,运行结果 x1 = -2.0074 y1 =2.2204e-016 x2 = -0.5198 y2 = 0 x3 =-0.5198 y3 =0 x4 =0.5993 y4 =-5.5511e-017 x5 = 1.6738 y5 =-1.1102e-016,（5）求在xx(3)附近的精确零点 从理论分析可知， 函数f(x)=(sin2x)e-0.1x-0.5|x| 0 是函数的一个零点。 但即便是以十分靠近该零点的值为搜索的初始值， 也找不到 ，而却找到了另一个零点。 原因是曲线没有穿越横轴。,三、求函数极值点,Matlab函数中只有处理极小值的指令。 求的是局部极值 1.一元函数极值 x,fval=fminbnd(fun,x1,x2) %x1、x2表示被研究区间的左右边界。 %x返回极小值点，fval返回该点的函数值。 2.多元函数极值 x,fval=fminsearch(fun,x0) %单纯形法求多元函数极值点指令，x0为初始值 x,fval=fminunc(fun,x0) %拟牛顿法求多元函数极值点指令，x0为初始值 %fun均是用字符串表示函数,x,fval=fminbnd(sin(x),-pi,pi) %求sinx在-,的极小值点 x = -1.5708 fval = -1.0000,x,fval=fminsearch( (x(1)-1)2+(x(2)-x(1)2)2,0,0) x = 1.0000 1.0000 fval = 2.6131e-010,x,fval=fminunc(x(1)2+x(2)2-x(1)*x(2)+2*x(1)-4*x(2),2,2) x = 0 2 fval = -4,四、数值积分,梯形方法trapz 辛普森公式quad 科特斯公式quad8 精度依次提高,梯形方法 matlab命令： x=eps:0.001:1; %eps表示什么？ y=sin(x)./x; %注意./ trapz(x,y) ans = 0.9461,由于辛普森公式quad、科特斯公式quad8这二种方法依 据的积分法不同于梯形法，因此它们的语法就和 trapz 不 同； 其语法为 quad(function,a,b) （quad8语法相同），其 中function是一已定义函数的名称（如sin, cos, sqrt, log 等或自己定义的函数），而 a, b是积分的下限和 上限。 和 trapz比较，quad, quad8不同之处在于这二者类似解 析式的积分式，只须设定上下限及定义要积分 的函数；而 trapz则是针对离散点型态的数据做积分。,quad(sin(x)./x,0,1) %注意./ Warning: Divide by zero. In inlineeval at 13 In inline.subsref at 25 In quad at 62 ans = 0.9461,quad(sin(x)./x,eps,1) ans = 0.9461,我们看一简单积分式 以下为 MATLAB 的程序 format long kq=quad(sqrt,0,1) kq = 0.66665956927202 kq8=quad8(sqrt,0,1) kq8 = 0.66666534581698,五、曲线的多项式拟合,1.多项式拟合 已知一些离散的数据点x=x0,x1,xm, 和对应的函数值y=y0,y1,ym，用n次多项式拟合这组数据 p,s=polyfit(x,y,n) %Fit polynomial to data. %拟合的n次多项式为p(x)=anxn+ an-1x n-1+a1x+a0 p返回的是一个向量p=an ,an-1,a1,a0 s是一个构架数组，其中s.normr表示误差,例：已知一组实验数据如下，求它的拟合曲线 首先画图观察 x=1 2 3 4 5 ; y=4 4.5 6 8 8.5; plot(x,y, r*) 点的形状比较接近一条直线，所以用一次多项式拟合,p,s=polyfit(x,y,1) p = 1.2500 2.4500 s = R: 2x2 double df: 3 normr: 0.8216 %下面画拟合曲线的图，首先取点 xn=1:0.001:5; yn=polyval(p,xn); %Evaluate polynomial. %若p=an ,an-1,a1,a0， %则此命令求多项式anxn+ an-1x n-1+a1x+a0在xn各分量处的函数值 plot(xn,yn,b,x,y,r*),2.其他形式,（1）已知一些离散的数据点x=x0,x1,xm,和对应的函数值y=y0,y1,ym，拟合成型如y=aebx 做变换lny=bx+lna,这是个一次多项式，可以用xi，lnyi拟合（i=0,1,m） 拟合的曲线为：t=bx+lna 用来拟合的数据：x=x0,x1,xm, t =lny0,lny1,lnym， 命令： x=x0,x1,xm; y=y0,y1,ym; t=log(y); p,s=polyfit(x,t,1); %此时p返回的是b lna 思考：如果想画y=aebx的图像该怎么做？,此时p返回的是b lna 思考：如果想画y=aebx的图像该怎么做？ 答案：xn=x0:步长:xm; %取一些点 %由于拟合的曲线是y=aebx %计算在xn各分量处的函数值 yn=exp(p(2)*exp(p(1)*xn); plot(x,y, r*,xn,yn, b) %画图 legend(数据点,拟合曲线) %最好加上图例标注,例：matlab命令： x=1.00,1.25,1.50,1.75,2.00; y=5.10,5.79,6.53,7.45,8.46; plot(x,y, r*) %通过观察，确定拟合曲线形如y=a*exp(b*x),离散数据点的图像，通过观察确定拟合曲线形如y=a*exp(b*x),t=log(y); p,s=polyfit(x,t,1) %p(1)=b,p(2)=lna,下面画图y=a*exp(b*x) xn=1:0.001:2; yn= exp(p(2)*exp(p(1)*xn); hold on plot(xn,yn) legend(数据点,拟合曲线) hold off,运行结果： p = 0.5057