数学实验课后习题解答

1、实验一曲线绘图【练习与思考】画出下列常见曲线的图形。以直角坐标方程表示的曲线：1. 立方曲线clear;x=-2:0.1:2;y=x.3;plot(x,y) 2. 立方抛物线clear;y=-2:0.1:2;x=y.3;plot(x,y)grid on 3. 高斯曲线clear;x=-3:0.1:3;y=exp(-x.2);plot(x,y);grid on%axis equal 以参数方程表示的曲线4. 奈尔抛物线clear;t=-3:0.05:3;x=t.3;y=t.2;plot(x,y)axis equalgrid on 5. 半立方抛物线clear;t=-3:0.05:3;x=t.2;

2、y=t.3;plot(x,y)%axis equalgrid on 6. 迪卡尔曲线clear;a=3;t=-6:0.1:6;x=3*a*t./(1+t.2);y=3*a*t.2./(1+t.2);plot(x,y) 7. 蔓叶线clear;a=3;t=-6:0.1:6;x=3*a*t.2./(1+t.2);y=3*a*t.3./(1+t.2);plot(x,y) 8. 摆线clear;clc;a=1;b=1;t=0:pi/50:6*pi;x=a*(t-sin(t);y=b*(1-cos(t);plot(x,y);axis equalgrid on 9. 内摆线（星形线）clear;a=1;t

3、=0:pi/50:2*pi;x=a*cos(t).3;y=a*sin(t).3;plot(x,y) 10. 圆的渐伸线（渐开线）clear;a=1;t=0:pi/50:6*pi;x=a*(cos(t)+t.*sin(t);y=a*(sin(t)+t.*cos(t);plot(x,y)grid on 11. 空间螺线cleara=3;b=2;c=1;t=0:pi/50:6*pi;x=a*cos(t);y=b*sin(t);z=c*t;plot3(x,y,z)grid on 以极坐标方程表示的曲线：12. 阿基米德线clear;a=1;phy=0:pi/50:6*pi;rho=a*phy;pola

4、r(phy,rho,'r-*') 13. 对数螺线clear;a=0.1;phy=0:pi/50:6*pi;rho=exp(a*phy);polar(phy,rho) 14. 双纽线clear;a=1;phy=-pi/4:pi/50:pi/4;rho=a*sqrt(cos(2*phy);polar(phy,rho)hold onpolar(phy,-rho) 15. 双纽线clear;a=1;phy=0:pi/50:pi/2;rho=a*sqrt(sin(2*phy);polar(phy,rho)hold onpolar(phy,-rho)16. 四叶玫瑰线clear;clos

5、ea=1;phy=0:pi/50:2*pi;rho=a*sin(2*phy);polar(phy,rho)17. 三叶玫瑰线clear;closea=1;phy=0:pi/50:2*pi;rho=a*sin(3*phy);polar(phy,rho) 18. 三叶玫瑰线clear;closea=1;phy=0:pi/50:2*pi;rho=a*cos(3*phy);polar(phy,rho) 实验二极限与导数【练习与思考】1 求下列各极限（1） （2） （3）clear;syms ny1=limit(1-1/n)n,n,inf)y2=limit(n3+3n)(1/n),n,inf)y3=li

6、mit(sqrt(n+2)-2*sqrt(n+1)+sqrt(n),n,inf) y1 =1/exp(1)y2 =3y3 =0 （4） （5） （6）clear;syms x ;y4=limit(2/(x2-1)-1/(x-1),x,1)y5=limit(x*cot(2*x),x,0)y6=limit(sqrt(x2+3*x)-x,x,inf) y4 =-1/2y5 =1/2y6 =3/2 （7） （8） （9）clear;syms x my7=limit(cos(m/x),x,inf)y8=limit(1/x-1/(exp(x)-1),x,1)y9=limit(1+x)(1/3)-1)/x,

7、x,0) y7 =1y8 =(exp(1) - 2)/(exp(1) - 1)y9 =1/3 2 考虑函数作出图形，并说出大致单调区间；使用diff求，并求确切的单调区间。clear;close;syms x;f=3*x2*sin(x3);ezplot(f,-2,2)grid on 大致的单调增区间：-2,-1.7,-1.3,1.2,1.7,2;大致的单点减区间：-1.7,-1.3,1.2,1.7 f1=diff(f,x,1)ezplot(f1,-2,2)line(-5,5,0,0)grid onaxis(-2.1,2.1,-60,120)f1 =6*x*sin(x3) + 9*x4*cos(

8、x3)用fzero函数找的零点，即原函数的驻点x1=fzero('6*x*sin(x3) + 9*x4*cos(x3)',-2,-1.7)x2=fzero('6*x*sin(x3) + 9*x4*cos(x3)',-1.7,-1.5)x3=fzero('6*x*sin(x3) + 9*x4*cos(x3)',-1.5,-1.1)x4=fzero('6*x*sin(x3) + 9*x4*cos(x3)',0)x5=fzero('6*x*sin(x3) + 9*x4*cos(x3)',1,1.5)x6=fzero(&#

9、39;6*x*sin(x3) + 9*x4*cos(x3)',1.5,1.7)x7=fzero('6*x*sin(x3) + 9*x4*cos(x3)',1.7,2)x1 = -1.9948x2 = -1.6926x3 = -1.2401x4 = 0x5 = 1.2401x6 = 1.6926x7 = 1.9948 确切的单调增区间：-1.9948,-1.6926,-1.2401,1.2401,1.6926,1.9948确切的单调减区间：-2,-1.9948,-1.6926,-1.2401,1.2401,1.6926,1.9948,23 对于下列函数完成下列工作，并写出

10、总结报告，评论极值与导数的关系，(i) 作出图形，观测所有的局部极大、局部极小和全局最大、全局最小值点的粗略位置；(iI) 求所有零点（即的驻点）；(iii) 求出驻点处的二阶导数值;(iv) 用fmin求各极值点的确切位置;(v) 局部极值点与有何关系?(1) (2) (3) clear;close;syms x;f=x2*sin(x2-x-2)ezplot(f,-2,2)grid on f =x2*sin(x2 - x - 2) 局部极大值点为：-1.6,局部极小值点为为：-0.75，-1.6全局最大值点为为：-1.6，全局最小值点为：-3f1=diff(f,x,1)ezplot(f1,-

11、2,2)line(-5,5,0,0)grid onaxis(-2.1,2.1,-6,20) f1 =2*x*sin(x2 - x - 2) + x2*cos(x2 - x - 2)*(2*x - 1)用fzero函数找的零点，即原函数的驻点x1=fzero('2*x*sin(x2-x-2)+x2*cos(x2-x-2)*(2*x-1)',-2,-1.2)x2=fzero('2*x*sin(x2-x-2)+x2*cos(x2-x-2)*(2*x-1)',-1.2,-0.5)x3=fzero('2*x*sin(x2-x-2)+x2*cos(x2-x-2)*(

12、2*x-1)',-0.5,1.2)x4=fzero('2*x*sin(x2-x-2)+x2*cos(x2-x-2)*(2*x-1)',1.2,2)x1 = -1.5326x2 = -0.7315x3 = -3.2754e-027x4 = 1.5951 ff=(x) x.2.*sin(x.2-x-2)ff(-2),ff(x1),ff(x2),ff(x3),ff(x4),ff(2) ff = (x)x.2.*sin(x.2-x-2)ans = -3.0272ans = 2.2364ans = -0.3582ans = -9.7549e-054ans = -2.2080ans

13、 = 0 实验三级数【练习与思考】1. 用taylor命令观测函数的Maclaurin展开式的前几项, 然后在同一坐标系里作出函数和它的Taylor展开式的前几项构成的多项式函数的图形，观测这些多项式函数的图形向的图形的逼近的情况(1) clear;syms xy=asin(x);y1=taylor(y,0,1)y2=taylor(y,0,5)y3=taylor(y,0,10)y4=taylor(y,0,15)x=-1:0.1:1;y=subs(y,x);y1=subs(y1,x);y2=subs(y2,x);y3=subs(y3,x);y4=subs(y4,x);plot(x,y,x,y1,

14、':',x,y2,'-.',x,y3,'-',x,y4,':','linewidth',3) y1 =0y2 =x3/6 + xy3 =(35*x9)/1152 + (5*x7)/112 + (3*x5)/40 + x3/6 + xy4 =(231*x13)/13312 + (63*x11)/2816 + (35*x9)/1152 + (5*x7)/112 + (3*x5)/40 + x3/6 + x (2) clear;syms xy=atan(x);y1=taylor(y,0,3)y2=taylor(y,0,5

15、),y3=taylor(y,0,10),y4=taylor(y,0,15)x=-1:0.1:1;y=subs(y,x);y1=subs(y1,x);y2=subs(y2,x);y3=subs(y3,x);y4=subs(y4,x);plot(x,y,x,y1,':',x,y2,'-.',x,y3,'-',x,y4,':','linewidth',3) y1 =xy2 =x - x3/3y3 =x9/9 - x7/7 + x5/5 - x3/3 + xy4 =x13/13 - x11/11 + x9/9 - x7/7

16、 + x5/5 - x3/3 + x(3) clear;syms xy=exp(x2);y1=taylor(y,0,3)y2=taylor(y,0,5)y3=taylor(y,0,10)y4=taylor(y,0,15)x=-1:0.1:1;y=subs(y,x);y1=subs(y1,x);y2=subs(y2,x);y3=subs(y3,x);y4=subs(y4,x);plot(x,y,x,y1,':',x,y2,'-.',x,y3,'-',x,y4,':','linewidth',3) y1 =x2 + 1

17、y2 =x4/2 + x2 + 1y3 =x8/24 + x6/6 + x4/2 + x2 + 1y4 =x14/5040 + x12/720 + x10/120 + x8/24 + x6/6 + x4/2 + x2 + 1(4) clear;syms xy=sin(x)2;y1=taylor(y,0,1)y2=taylor(y,0,5)y3=taylor(y,0,10)y4=taylor(y,0,15)x=-pi:0.1:pi;y=subs(y,x);y1=subs(y1,x);y2=subs(y2,x);y3=subs(y3,x);y4=subs(y4,x);plot(x,y,x,y1,&

18、#39;:',x,y2,'-.',x,y3,'-',x,y4,':','linewidth',3) y1 =0y2 =x2 - x4/3y3 =- x8/315 + (2*x6)/45 - x4/3 + x2y4 =(4*x14)/42567525 - (2*x12)/467775 + (2*x10)/14175 - x8/315 + (2*x6)/45 - x4/3 + x2 (5) clear;syms xy=exp(x)/(1-x);y1=taylor(y,0,3)y2=taylor(y,0,5)y3=taylor(

19、y,0,10)y4=taylor(y,0,15)x=-1:0.1:0;y=subs(y,x);y1=subs(y1,x);y2=subs(y2,x);y3=subs(y3,x);y4=subs(y4,x);plot(x,y,x,y1,':',x,y2,'-.',x,y3,'-',x,y4,':','linewidth',3) y1 =(5*x2)/2 + 2*x + 1y2 =(65*x4)/24 + (8*x3)/3 + (5*x2)/2 + 2*x + 1y3 =(98641*x9)/36288 + (1096

20、01*x8)/40320 + (685*x7)/252 + (1957*x6)/720 + (163*x5)/60 + (65*x4)/24 + (8*x3)/3 + (5*x2)/2 + 2*x + 1y4 =(47395032961*x14)+ (8463398743*x13)/3113510400 + (260412269*x12)/95800320 + (13563139*x11)/4989600 + (9864101*x10)/3628800 + (98641*x9)/36288 + (109601*x8)/40320 + (685*x7)/252 + (1

21、957*x6)/720 + (163*x5)/60 + (65*x4)/24 + (8*x3)/3 + (5*x2)/2 + 2*x + 1 (6) clear;syms xy=log(x+sqrt(1+x2);y1=taylor(y,0,3)y2=taylor(y,0,5)y3=taylor(y,0,10)y4=taylor(y,0,15)x=-1:0.1:1;y=subs(y,x);y1=subs(y1,x);y2=subs(y2,x);y3=subs(y3,x);y4=subs(y4,x);plot(x,y,x,y1,':',x,y2,'-.',x,y3,

22、'-',x,y4,':','linewidth',3) y1 =xy2 =x - x3/6y3 =(35*x9)/1152 - (5*x7)/112 + (3*x5)/40 - x3/6 + xy4 =(231*x13)/13312 - (63*x11)/2816 + (35*x9)/1152 - (5*x7)/112 + (3*x5)/40 - x3/6 + x 2. 求公式中的数的值.k=4 5 6 7 8;syms nsymsum(1./n.(2*k),1,inf) ans = pi8/9450, pi10/93555, (691*pi12

23、)/638512875, (2*pi14)/18243225, (3617*pi16)/325641566250 3. 利用公式来计算的近似值。精确到小数点后100位,这时应计算到这个无穷级数的前多少项?请说明你的理由.解：Matlab代码为clear;clc;closeepsl=1.0e-100;ep=1;fn=1;a=1;n=1;while ep>epsla=a+fn;n=n+1;fn=fn/n;ep=fn;endfnvpa(a,100)n fn = 8.3482e-101ans =2.7182818284590455348848081484902650117874145507812

24、5n = 70 精确到小数点后100位,这时应计算到这个无穷级数的前71项,理由是误差小于10的负100次方，需要最后一项小于10的负100次方，由上述循环知n=70时最后一项小于10的负100次方，故应计算到这个无穷级数的前71项.4. 用练习3中所用观测法判断下列级数的敛散性(1) clear;clc;epsl=0.000001;N=50000;p=1000;syms nUn=1/(n2+n3);s1=symsum(Un,1,N);s2=symsum(Un,1,N+p);sa=vpa(s2-s1);sa=setstr(sa);sa=str2num(sa);fprintf('级数&#

25、39;)disp(Un)if sa<epsl disp('收敛')else disp('发散')end 级数1/(n3 + n2)收敛 clear;closesyms ns=;for k=1:100s(k)=symsum(1/(n3 + n2),1,k);endplot(s,'.') (2) clear;clc;epsl=0.000001;N=50000;p=1000;syms nUn=1/(n*2n);s1=symsum(Un,1,N);s2=symsum(Un,1,N+p);sa=vpa(s2-s1);sa=setstr(sa);sa=

26、str2num(sa);fprintf('级数')disp(Un)if sa<epsl disp('收敛')else disp('发散')end 级数1/(2n*n)收敛 clear;closesyms ns=;for k=1:100s(k)=symsum(1/(2n*n),1,k);endplot(s,'.') (3) clear;clc;epsl=0.00000000000001;N=50000;p=100;syms nUn=1/sin(n);s1=symsum(Un,1,N);s2=symsum(Un,1,N+p);s

27、a=vpa(s2-s1);sa=setstr(sa);sa=str2num(sa);fprintf('级数')disp(Un)if abs(sa)<epsl disp('收敛')else disp('发散')end 级数1/sin(n)发散 clear;closesyms ns=;for k=1:100s(k)=symsum(1/sin(n),1,k);endplot(s,'.') 发散 (4) clear;clc;epsl=0.0000001;N=50000;p=1000;syms nUn=log(n)/(n3);s1=s

28、ymsum(Un,1,N);s2=symsum(Un,1,N+p);sa=vpa(s2-s1);sa=setstr(sa);sa=str2num(sa);fprintf('级数')disp(Un)if sa<epsl disp('收敛')else disp('发散')end 级数log(n)/n3收敛 clear;closesyms ns=;for k=1:100s(k)=symsum(log(n)/n3,1,k);endplot(s,'.') (5) clear;closesyms ns=;he=0;for k=1:100

29、he=he+factorial(k)/kk;s(k)=he;endplot(s,'.') (6) clear;clc;epsl=0.0000001;N=50000;p=1000;syms nUn=1/log(n)n;s1=symsum(Un,3,N);s2=symsum(Un,3,N+p);sa=vpa(s2-s1);sa=setstr(sa);sa=str2num(sa);fprintf('级数')disp(Un)if sa<epsl disp('收敛')else disp('发散')end 级数1/log(n)n收敛 c

30、lear;closesyms ns=;for k=3:100s(k)=symsum(1/log(n)n,3,k);endplot(s,'.') (7) clear;clc;epsl=0.0000001;N=50000;p=100;syms nUn=1/(log(n)*n);s1=symsum(Un,3,N);s2=symsum(Un,3,N+p);sa=vpa(s2-s1);sa=setstr(sa);sa=str2num(sa);fprintf('级数')disp(Un)if (sa)<epsl disp('收敛')else disp(&

31、#39;发散')end 级数1/(n*log(n)发散 clear;closesyms ns=;for k=3:300s(k)=symsum(1/(n*log(n),2,k);endplot(s,'.') (8) clear;clc;epsl=0.0000001;N=50000;p=100;syms nUn=(-1)n*n/(n2+1);s1=symsum(Un,3,N);s2=symsum(Un,3,N+p);sa=vpa(s2-s1);sa=setstr(sa);sa=str2num(sa);fprintf('级数')disp(Un)if (sa)&

32、lt;epsl disp('收敛')else disp('发散')end 级数(-1)n*n)/(n2 + 1)收敛 clear;closesyms ns=;for k=3:300s(k)=symsum(-1)n*n/(n2+1),2,k);endplot(s,'.') 实验四积分【练习与思考】1(不定积分)用int计算下列不定积分，并用diff验证，解：Matlab代码为：syms xy1=x*sin(x2);y2=1/(1+cos(x);y3=1/(exp(x)+1);y4=asin(x);y5=sec(x)3;f1=int(y1)f2=in

33、t(y2)f3=int(y3)f4=int(y4)f5=int(y5) dy=simplify(diff(f1;f2;f3;f4;f5) dy = x*sin(x2) tan(x/2)2/2 + 1/2 1/(exp(x) + 1) asin(x) (cot(pi/4 + x/2)*(tan(pi/4 + x/2)2/2 + 1/2)/2 + 1/(2*cos(x) + tan(x)2/cos(x) f1 =-cos(x2)/2f2 =tan(x/2)f3 =x - log(exp(x) + 1)f4 =x*asin(x) + (1 - x2)(1/2)f5 =log(tan(pi/4 + x

34、/2)/2 + tan(x)/(2*cos(x) 2(定积分)用trapz,quad,int计算下列定积分,解：Matlab代码为clear;x=(0+eps):0.05:1;y1=sin(x)./x;f1=trapz(x,y1) f1 =0.9460 fun1=(x)sin(x)./x;f12=quad(fun1,0+eps,1) f12 = 0.9461 f13=vpa(int('sin(x)/x',0,1),5) f13 =0.94608 3(椭圆的周长) 用定积分的方法计算椭圆的周长解：椭圆的参数方程为由参数曲线的弧长公式得Matlab代码为s=vpa(int('

35、;sqrt(5*sin(t)2+4)','t',0,2*pi),5) s =15.865 4（二重积分）计算数值积分解：fxy=(x,y)1+x+y;ylow=(x)1-sqrt(1-x.2);yup=(x)1+sqrt(1-x.2);s=quad2d(fxy,-1,1,ylow,yup) s =6.2832 或符号积分法：syms x yxi=int(1+x+y,y,1-sqrt(1-x2),1+sqrt(1-x2);s=int(xi,x,-1,1) s =2*pi 5（假奇异积分）用trapz,quad8计算积分，会出现什么问题？分析原因，并求出正确的解。解：Mat

36、lab代码为clearx=-1:0.05:1;y=x.(1/3).*cos(x);s1=trapz(x,y)fun5=(x)x.(1/3).*cos(x);s2=quad(fun5,-1,1)int('x(1/3)*cos(x)','x',-1,1) s1 = 0.9036 + 0.5217is2 = 0.9114 + 0.5262iWarning: Explicit integral could not be found. ans =int(x(1/3)*cos(x), x = -1.1) ，原函数不存在，不能用int函数运算。用梯形法和辛普森法计算数值积分时

37、，由于对负数的开三次方运算结果为复数，所以导致结果错误且为复数；显然被积函数为奇函数，在对称区间上的积分等于0，此时可以这样处理：（1）重新定义被积函数%fun5.mfunction y=fun5(x)m,n=size(x);for k=1:mfor l=1:ny(k,l)=nthroot(x(k,l),3)*cos(x(k,l);endendend用辛普森法：s=quad('fun5',-1,1) s = 0 用梯形法clear;x=-1:0.01:1;y=fun5(x);s=trapz(x,y) s = -1.3878e-017 6（假收敛现象）考虑积分，（1）用解析法求;

38、clear;syms x k;Ik=int(abs(sin(x),0,k*pi) Warning: Explicit integral could not be found. Ik =int(abs(sin(x), x = 0.pi*k) （2）分别用trapz,quad和quad8求和,发现什么问题？clear;for k=4:2:8;x=0:pi/1000:k*pi;y=abs(sin(x);trapz(x,y)end ans = 8.0000ans = 12.0000ans = 16.0000 for k=4:2:8fun6=(x)abs(sin(x);quad(fun6,0,k*pi)

39、end ans = 8.0000ans = 12.0000ans = 16.0000 7(Simpson积分法)编制一个定步长Simpson法数值积分程序.计算公式为其中为偶数,解：Matlab代码为%fun7.mfunction y=fun7(f_name,a,b,n)%f_name为被积函数%a,b为积分区间%n为偶数，用来确定步长h=(b-a)/nif mod(n,2)=0 disp('n必须为偶数') return;endif nargin<4n=100;endif nargin<3disp('请输入积分区间')endif nargin=0d

40、isp('error')endh=(b-a)/n;x=a:h:b;s=0;for k=1:n+1if k=1|k=(n+1) xishu=1;elseif mod(k,2)=0 xishu=4;else xishu=2;ends=s+feval(f_name,x(k)*xishu;endy=s*h/3;end8(广义积分)计算广义积分,并验证公式.解：Matlab代码为clear;syms xs1=vpa(int(exp(-x2)/(1+x4),-inf,inf),5)s2=quad(x)tan(x)./sqrt(x),0,1)s3=quad(x)sin(x)./sqrt(1-

41、x.2),0,1)s4=vpa(int(exp(-x2/2)/sqrt(2*pi),-inf,inf),5)s5=int(sin(x)./x,0+eps,inf) s1 =1.4348s2 = 0.7968s3 = 0.8933s4 =1.0s5 =pi/2 - sinint(1/4503599627370496) 实验五二元函数的图形【练习与思考】1. 画出空间曲线在范围内的图形，并画出相应的等高线。clear;x=-30:0.5:30;y=-30:0.5:30;X,Y=meshgrid(x,y);Z=10*sin(sqrt(X.2+Y.2)./sqrt(1+X.2+Y.2);mesh(X,

42、Y,Z) 2. 根据给定的参数方程，绘制下列曲面的图形。a) 椭球面clear;u=0:pi/50:2*pi;v=0:pi/50:pi;U,V=meshgrid(u,v);x=3*cos(U).*sin(V);y=2*cos(U).*cos(V);z=sin(U);mesh(x,y,z) b) 椭圆抛物面clear;u=0:pi/50:pi/4;v=0:pi/50:2*pi;U,V=meshgrid(u,v);x=3*U.*sin(V);y=2*U.*cos(V);z=4*U.2;mesh(x,y,z)axis equal c) 单叶双曲面clear;u=0:pi/15:pi;v=0:pi/1

43、5:2*pi;U,V=meshgrid(u,v);x=3*sec(U).*sin(V);y=2*sec(U).*cos(V);z=4*tan(U);mesh(x,y,z) d) 双曲抛物面clearu=-3:0.1:3;U,V=meshgrid(u);x=U;y=V;z=(U.2-V.2)/3;mesh(x,y,z) e) 旋转面clear;u=1:0.1:5;v=0:pi/30:2*pi;U,V=meshgrid(u,v);x=log(U).*sin(V);y=log(U).*cos(V);z=U;mesh(x,y,z)axis equal f) 圆锥面clear;u=-5:0.1:5;v=

44、0:pi/30:2*pi;U,V=meshgrid(u,v);x=(U).*sin(V);y=(U).*cos(V);z=U;mesh(x,y,z)axis equal g) 环面clear;u=0:pi/30:2*pi;v=u;U,V=meshgrid(u,v);x=(3+0.4*cos(U).*cos(V);y=(3+0.4*cos(U).*sin(V);z=0.4*sin(V);mesh(x,y,z) h) 正螺面clear;u=0:pi/30:pi;v=0:pi/30:10*pi;U,V=meshgrid(u,v);x=U.*sin(V);y=U.*cos(V);z=4*V;mesh(

45、x,y,z)colorbar 3. 在一丘陵地带测量搞程,x和y方向每隔100米测一个点,得高程见表5-2,试拟合一曲面,确定合适的模型,并由此找出最高点和该点的高程.表5-2 高程数据y x100200300400100200300400636698680662697712674626624630598552478478412334clc;clear;x1=100 100 100 100 200 200 200 200 300 300 300 300 400 400 400 400;x2=100 200 300 400 100 200 300 400 100 200 300 400 100

46、200 300 400;y=636 698 680 662 697 712 674 626 624 630 598 552 478 478 412 334'x=x1',x2'x0=1 1 1 1 1;beta=lsqcurvefit('heigh',x0,x,y)%绘图： a1=100:5:400;a2=a1;xx1,xx2=meshgrid(a1,a2);Z=beta(1)+beta(2)*xx1+beta(3)*xx2+beta(4)*xx1.2+beta(5)*xx2.2;mesh(xx1,xx2,Z)Local minimum possible.

47、lsqcurvefit stopped because the final change in the sum of squares relative to its initial value is less than the default value of the function tolerance.beta = Columns 1 through 5 538.4375 1.4901 0.6189 -0.0046 -0.0017 contour(xx1,xx2,Z,30),colorbar%计算最高点及高程x0=100,100;options=optimset('largesca

48、le','off');%设置下界lb=0,0;%无上界ub=;x,fval=fmincon('height',x0,lb,ub,options) Warning: Options LargeScale = 'off' and Algorithm = 'trust-region-reflective' conflict.Ignoring Algorithm and running active-set algorithm. To run trust-region-reflective, setLargeScale = 

49、9;on'. To run active-set without this warning, use Algorithm = 'active-set'.> In fmincon at 445Local minimum possible. Constraints satisfied.fmincon stopped because the predicted change in the objective functionis less than the default value of the function tolerance and constraints w

50、ere satisfied to within the default value of the constraint tolerance.No active inequalities.x = 161.9676 182.0320fval = -715.4403 heigh和height两个函数分别定义如下：(应写在m文件中)%heigh.mfunction f=heigh(beta,xdata)xx1=xdata(:,1);xx2=xdata(:,2);f=beta(1)+beta(2)*xx1+beta(3)*xx2+beta(4)*xx1.2+beta(5)*xx2.*xx1+beta(6

51、)*xx2.2;end%height.mfunction y=height(x)y=-(434.0000+1.9079*x(1)+1.0366*x(2)-0.0017*x(1).2-0.0046*x(2).*x(1)-0.0017*x(2).2);end实验六多元函数的极值【练习与思考】1.求的极值，并对图形进行观测。解：Maltab代码为syms x y;z=x4+y4-4*x*y+1;dzx=diff(z,x);dzy=diff(z,y);x,y=solve(dzx,dzy,x,y) x = 0 1 -1 (-1)(3/4) -(-1)(3/4) -i i -(-1)(3/4)*i (-1)(3/4)*iy = 0 1 -1 (-1)(1/4) -(-1)(1/4) i -i (-1)(1/4)*i -(-1)(1/4)*i 经计算可知，函数的驻点为(0,0)、(1,1)、(-1,-1)ezmeshc(z,-2