MATLAB数学实验第二版答案 胡良剑
(完整版)MATLAB数学实验第二版答案(胡良剑)
数学实验答案Chapter 1Page20,ex1(5) 等于[exp(1),exp(2);exp(3),exp(4)](7) 3=1*3, 8=2*4(8) a为各列最小值,b为最小值所在的行号(10) 1>=4,false, 2>=3,false, 3>=2, ture, 4>=1,ture(11) 答案表明:编址第2元素满足不等式(30>=20)和编址第4元素满足不等式(40>=10)(12) 答案表明:编址第2行第1列元素满足不等式(30>=20)和编址第2行第2列元素满足不等式(40>=10)Page20, ex2(1)a, b, c的值尽管都是1,但数据类型分别为数值,字符,逻辑,注意a与c相等,但他们不等于b(2)double(fun)输出的分别是字符a,b,s,(,x,)的ASCII码Page20,ex3>> r=2;p=0.5;n=12;>> T=log(r)/n/log(1+0.01*p)Page20,ex4>> x=-2:0.05:2;f=x.^4-2.^x;>> [fmin,min_index]=min(f)最小值最小值点编址>> x(min_index)ans =0.6500 最小值点>> [f1,x1_index]=min(abs(f)) 求近似根--绝对值最小的点f1 =0.0328x1_index =24>> x(x1_index)ans =-0.8500>> x(x1_index)=[];f=x.^4-2.^x; 删去绝对值最小的点以求函数绝对值次小的点>> [f2,x2_index]=min(abs(f)) 求另一近似根--函数绝对值次小的点f2 =0.0630x2_index =65>> x(x2_index)ans =1.2500>> z=magic(10)z =92 99 1 8 15 67 74 51 58 4098 80 7 14 16 73 55 57 64 414 81 88 20 22 54 56 63 70 4785 87 19 21 3 60 62 69 71 2886 93 25 2 9 61 68 75 52 3417 24 76 83 90 42 49 26 33 6523 5 82 89 91 48 30 32 39 6679 6 13 95 97 29 31 38 45 7210 12 94 96 78 35 37 44 46 5311 18 100 77 84 36 43 50 27 59>> sum(z)>> sum(diag(z))>> z(:,2)/sqrt(3)>> z(8,:)=z(8,:)+z(3,:)Chapter 2Page 45 ex1先在编辑器窗口写下列M函数,保存为eg2_1.m function [xbar,s]=ex2_1(x)n=length(x);xbar=sum(x)/n;s=sqrt((sum(x.^2)-n*xbar^2)/(n-1));例如>>x=[81 70 65 51 76 66 90 87 61 77];>>[xbar,s]=ex2_1(x)Page 45 ex2s=log(1);n=0;while s<=100n=n+1;s=s+log(1+n);endm=nPage 40 ex3clear;F(1)=1;F(2)=1;k=2;x=0;e=1e-8; a=(1+sqrt(5))/2;while abs(x-a)>ek=k+1;F(k)=F(k-1)+F(k-2); x=F(k)/F(k-1);enda,x,k计算至k=21可满足精度clear;tic;s=0;for i=1:1000000s=s+sqrt(3)/2^i;ends,toctic;s=0;i=1;while i<=1000000s=s+sqrt(3)/2^i;i=i+1;ends,toctic;s=0;i=1:1000000;s=sqrt(3)*sum(1./2.^i);s,tocPage 45 ex5t=0:24;c=[15 14 14 14 14 15 16 18 20 22 23 25 28 ...31 32 31 29 27 25 24 22 20 18 17 16];plot(t,c)Page 45 ex6(1)x=-2:0.1:2;y=x.^2.*sin(x.^2-x-2);plot(x,y)y=inline('x^2*sin(x^2-x-2)');fplot(y,[-2 2]) (2)参数方法t=linspace(0,2*pi,100);x=2*cos(t);y=3*sin(t); plot(x,y)(3)x=-3:0.1:3;y=x;[x,y]=meshgrid(x,y);z=x.^2+y.^2;surf(x,y,z)(4)x=-3:0.1:3;y=-3:0.1:13;[x,y]=meshgrid(x,y);z=x.^4+3*x.^2+y.^2-2*x-2*y-2*x.^2.*y+6;surf(x,y,z)(5)t=0:0.01:2*pi;x=sin(t);y=cos(t);z=cos(2*t);plot3(x,y,z)(6)theta=linspace(0,2*pi,50);fai=linspace(0,pi/2,20); [theta,fai]=meshgrid(theta,fai);x=2*sin(fai).*cos(theta);y=2*sin(fai).*sin(theta);z=2*cos(fai);surf(x,y,z)(7)x=linspace(0,pi,100);y1=sin(x);y2=sin(x).*sin(10*x);y3=-sin(x);plot(x,y1,x,y2,x,y3)page45, ex7x=-1.5:0.05:1.5;y=1.1*(x>1.1)+x.*(x<=1.1).*(x>=-1.1)-1.1*(x<-1.1);plot(x,y)page45,ex9clear;close;x=-2:0.1:2;y=x;[x,y]=meshgrid(x,y);a=0.5457;b=0.7575;p=a*exp(-0.75*y.^2-3.75*x.^2-1.5*x).*(x+y>1);p=p+b*exp(-y.^2-6*x.^2).*(x+y>-1).*(x+y<=1);p=p+a*exp(-0.75*y.^2-3.75*x.^2+1.5*x).*(x+y<=-1);mesh(x,y,p)page45, ex10lookfor lyapunovhelp lyap>> A=[1 2 3;4 5 6;7 8 0];C=[2 -5 -22;-5 -24 -56;-22 -56 -16];>> X=lyap(A,C)X =1.0000 -1.0000 -0.0000-1.0000 2.0000 1.0000-0.0000 1.0000 7.0000Chapter 3Page65 Ex1>> a=[1,2,3];b=[2,4,3];a./b,a.\b,a/b,a\bans =0.5000 0.5000 1.0000ans =2 2 1ans =0.6552 一元方程组x[2,4,3]=[1,2,3]的近似解ans =0 0 00 0 00.6667 1.3333 1.0000矩阵方程[1,2,3][x11,x12,x13;x21,x22,x23;x31,x32,x33]=[2,4,3]的特解Page65 Ex 2(1)>> A=[4 1 -1;3 2 -6;1 -5 3];b=[9;-2;1];>> rank(A), rank([A,b]) [A,b]为增广矩阵ans =3ans =3 可见方程组唯一解>> x=A\bx =2.38301.48942.0213(2)>> A=[4 -3 3;3 2 -6;1 -5 3];b=[-1;-2;1];>> rank(A), rank([A,b])ans =3ans =3 可见方程组唯一解>> x=A\bx =-0.4706-0.2941(3)>> A=[4 1;3 2;1 -5];b=[1;1;1];>> rank(A), rank([A,b])ans =2ans =3 可见方程组无解>> x=A\bx =0.3311-0.1219 最小二乘近似解(4)>> a=[2,1,-1,1;1,2,1,-1;1,1,2,1];b=[1 2 3]';%注意b的写法>> rank(a),rank([a,b])ans =3ans =3 rank(a)==rank([a,b])<4说明有无穷多解>> a\bans =110 一个特解Page65 Ex3>> a=[2,1,-1,1;1,2,1,-1;1,1,2,1];b=[1,2,3]';>> x=null(a),x0=a\bx =-0.62550.6255-0.20850.4170x0 =11通解kx+x0Page65 Ex 4>> x0=[0.2 0.8]';a=[0.99 0.05;0.01 0.95];>> x1=a*x, x2=a^2*x, x10=a^10*x>> x=x0;for i=1:1000,x=a*x;end,xx =0.83330.1667>> x0=[0.8 0.2]';>> x=x0;for i=1:1000,x=a*x;end,xx =0.83330.1667>> [v,e]=eig(a)v =0.9806 -0.70710.1961 0.7071e =1.0000 00 0.9400>> v(:,1)./xans =1.17671.1767 成比例,说明x是最大特征值对应的特征向量Page65 Ex5用到公式(3.11)(3.12)>> B=[6,2,1;2.25,1,0.2;3,0.2,1.8];x=[25 5 20]'; >> C=B/diag(x)C =0.2400 0.4000 0.05000.0900 0.2000 0.01000.1200 0.0400 0.0900>> A=eye(3,3)-CA =0.7600 -0.4000 -0.0500-0.0900 0.8000 -0.0100-0.1200 -0.0400 0.9100>> D=[17 17 17]';x=A\Dx =37.569625.786224.7690Page65 Ex 6(1)>> a=[4 1 -1;3 2 -6;1 -5 3];det(a),inv(a),[v,d]=eig(a) ans =-94ans =0.2553 -0.0213 0.04260.1596 -0.1383 -0.22340.1809 -0.2234 -0.0532v =0.0185 -0.9009 -0.3066-0.7693 -0.1240 -0.7248-0.6386 -0.4158 0.6170d =-3.0527 0 00 3.6760 00 0 8.3766(2)>> a=[1 1 -1;0 2 -1;-1 2 0];det(a),inv(a),[v,d]=eig(a) ans =1ans =2.0000 -2.0000 1.00001.0000 -1.0000 1.00002.0000 -3.0000 2.0000v =-0.5773 0.5774 + 0.0000i 0.5774 - 0.0000i-0.5773 0.5774 0.5774-0.5774 0.5773 - 0.0000i 0.5773 + 0.0000id =1.0000 0 00 1.0000 + 0.0000i 00 0 1.0000 - 0.0000i(3)>> A=[5 7 6 5;7 10 8 7;6 8 10 9;5 7 9 10]A =5 76 57 10 8 76 8 10 95 7 9 10>> det(A),inv(A), [v,d]=eig(A)ans =1ans =68.0000 -41.0000 -17.0000 10.0000-41.0000 25.0000 10.0000 -6.0000-17.0000 10.0000 5.0000 -3.000010.0000 -6.0000 -3.0000 2.0000v =0.8304 0.0933 0.3963 0.3803-0.5016 -0.3017 0.6149 0.5286-0.2086 0.7603 -0.2716 0.55200.1237 -0.5676 -0.6254 0.5209d =0.0102 0 0 00 0.8431 0 00 0 3.8581 00 0 0 30.2887(4)(以n=5为例)方法一(三个for)n=5;for i=1:n, a(i,i)=5;endfor i=1:(n-1),a(i,i+1)=6;endfor i=1:(n-1),a(i+1,i)=1;enda方法二(一个for)n=5;a=zeros(n,n);a(1,1:2)=[5 6];for i=2:(n-1),a(i,[i-1,i,i+1])=[1 5 6];enda(n,[n-1 n])=[1 5];a方法三(不用for)n=5;a=diag(5*ones(n,1));b=diag(6*ones(n-1,1));c=diag(ones(n-1,1));a=a+[zeros(n-1,1),b;zeros(1,n)]+[zeros(1,n);c,zeros(n-1,1)] 下列计算>> det(a)ans =665>> inv(a)ans =0.3173 -0.5865 1.0286 -1.6241 1.9489-0.0977 0.4887 -0.8571 1.3534 -1.62410.0286 -0.1429 0.5429 -0.8571 1.0286-0.0075 0.0376 -0.1429 0.4887 -0.58650.0015 -0.0075 0.0286 -0.0977 0.3173>> [v,d]=eig(a)v =-0.7843 -0.7843 -0.9237 0.9860 -0.92370.5546 -0.5546 -0.3771 -0.0000 0.3771-0.2614 -0.2614 0.0000 -0.1643 0.00000.0924 -0.0924 0.0628 -0.0000 -0.0628-0.0218 -0.0218 0.0257 0.0274 0.0257d =0.7574 0 0 0 00 9.2426 0 0 00 0 7.4495 0 00 0 0 5.0000 00 0 0 0 2.5505Page65 Ex 7(1)>> a=[4 1 -1;3 2 -6;1 -5 3];[v,d]=eig(a)v =0.0185 -0.9009 -0.3066-0.7693 -0.1240 -0.7248-0.6386 -0.4158 0.6170d =-3.0527 0 00 3.6760 00 0 8.3766>> det(v)ans =-0.9255 %v行列式正常, 特征向量线性相关,可对角化>> inv(v)*a*v 验算ans =-3.0527 0.0000 -0.00000.0000 3.6760 -0.0000-0.0000 -0.0000 8.3766>> [v2,d2]=jordan(a) 也可用jordanv2 =0.0798 0.0076 0.91270.1886 -0.3141 0.1256-0.1605 -0.2607 0.4213 特征向量不同d2 =8.3766 0 00 -3.0527 - 0.0000i 00 0 3.6760 + 0.0000i>> v2\a*v2ans =8.3766 0 0.00000.0000 -3.0527 0.00000.0000 0.0000 3.6760>> v(:,1)./v2(:,2) 对应相同特征值的特征向量成比例ans =2.44912.44912.4491(2)>> a=[1 1 -1;0 2 -1;-1 2 0];[v,d]=eig(a)v =-0.5773 0.5774 + 0.0000i 0.5774 - 0.0000i-0.5773 0.5774 0.5774-0.5774 0.5773 - 0.0000i 0.5773 + 0.0000id =1.0000 0 00 1.0000 + 0.0000i 00 0 1.0000 - 0.0000i>> det(v)ans =-5.0566e-028 -5.1918e-017i v的行列式接近0, 特征向量线性相关,不可对角化>> [v,d]=jordan(a)v =1 0 11 0 01 -1 0d =1 1 00 1 10 0 1 jordan标准形不是对角的,所以不可对角化(3)>> A=[5 7 6 5;7 10 8 7;6 8 10 9;5 7 9 10]A =5 76 57 10 8 76 8 10 95 7 9 10>> [v,d]=eig(A)v =0.8304 0.0933 0.3963 0.3803-0.5016 -0.3017 0.6149 0.5286-0.2086 0.7603 -0.2716 0.55200.1237 -0.5676 -0.6254 0.5209d =0.0102 0 0 00 0.8431 0 00 0 3.8581 00 0 0 30.2887>> inv(v)*A*vans =0.0102 0.0000 -0.0000 0.00000.0000 0.8431 -0.0000 -0.0000-0.0000 0.0000 3.8581 -0.0000-0.0000 -0.0000 0 30.2887本题用jordan不行, 原因未知(4)参考6(4)和7(1)Page65 Exercise 8只有(3)对称, 且特征值全部大于零, 所以是正定矩阵. Page65 Exercise 9(1)>> a=[4 -3 1 3;2 -1 3 5;1 -1 -1 -1;3 -2 3 4;7 -6 -7 0]>> rank(a)ans =3>> rank(a(1:3,:))ans =2>> rank(a([1 2 4],:)) 1,2,4行为最大无关组ans =3>> b=a([1 2 4],:)';c=a([3 5],:)';>> b\c 线性表示的系数ans =0.5000 5.0000-0.5000 1.00000 -5.0000Page65 Exercise 10>> a=[1 -2 2;-2 -2 4;2 4 -2]>> [v,d]=eig(a)v =0.3333 0.9339 -0.12930.6667 -0.3304 -0.6681-0.6667 0.1365 -0.7327d =-7.0000 0 00 2.0000 00 0 2.0000>> v'*vans =1.0000 0.0000 0.00000.0000 1.0000 00.0000 0 1.0000 v确实是正交矩阵Page65 Exercise 11设经过6个电阻的电流分别为i1, ..., i6. 列方程组如下20-2i1=a; 5-3i2=c; a-3i3=c; a-4i4=b; c-5i5=b; b-3i6=0;i1=i3+i4;i5=i2+i3;i6=i4+i5;计算如下>> A=[1 0 0 2 0 0 0 0 0;0 0 1 0 3 0 0 0 0;1 0 -1 0 0 -3 0 0 0; 1 -1 0 0 0 0 -4 0 0;0 -1 1 0 0 0 0 -5 0;0 1 0 0 0 0 0 0 -3; 0 0 0 1 0 -1 -1 0 0;0 0 0 0 -1 -1 0 1 0;0 0 0 0 0 0 -1 -1 1];>>b=[20 5 0 0 0 0 0 0 0]'; A\bans =13.34536.44018.54203.3274-1.18071.60111.72630.42042.1467Page65 Exercise 12>> A=[1 2 3;4 5 6;7 8 0];>> left=sum(eig(A)), right=sum(trace(A))left =6.0000right =6>> left=prod(eig(A)), right=det(A) 原题有错, (-1)^n应删去left =27.0000right =27>> fA=(A-p(1)*eye(3,3))*(A-p(2)*eye(3,3))*(A-p(3)*eye(3,3)) fA =1.0e-012 *0.0853 0.1421 0.02840.1421 0.1421 0-0.0568 -0.1137 0.1705>> norm(fA) f(A)范数接近0ans =2.9536e-013Chapter 4Page84 Exercise 1(1)roots([1 1 1])(2)roots([3 0 -4 0 2 -1])(3)p=zeros(1,24);p([1 17 18 22])=[5 -6 8 -5];roots(p)(4)p1=[2 3];p2=conv(p1, p1);p3=conv(p1, p2);p3(end)=p3(end)-4; %原p3最后一个分量-4roots(p3)Page84 Exercise 2fun=inline('x*log(sqrt(x^2-1)+x)-sqrt(x^2-1)-0.5*x');fzero(fun,2)Page84 Exercise 3fun=inline('x^4-2^x');fplot(fun,[-2 2]);grid on;fzero(fun,-1),fzero(fun,1),fminbnd(fun,0.5,1.5)Page84 Exercise 4fun=inline('x*sin(1/x)','x');fplot(fun, [-0.1 0.1]);x=zeros(1,10);for i=1:10, x(i)=fzero(fun,(i-0.5)*0.01);end;x=[x,-x]Page84 Exercise 5fun=inline('[9*x(1)^2+36*x(2)^2+4*x(3)^2-36;x(1)^2-2*x(2)^2-20*x(3);16*x(1)-x(1)^3-2*x(2)^ 2-16*x(3)^2]','x');[a,b,c]=fsolve(fun,[0 0 0])Page84 Exercise 6fun=@(x)[x(1)-0.7*sin(x(1))-0.2*cos(x(2)),x(2)-0.7*cos(x(1))+0.2*sin(x(2))];[a,b,c]=fsolve(fun,[0.5 0.5])Page84 Exercise 7clear; close; t=0:pi/100:2*pi;x1=2+sqrt(5)*cos(t); y1=3-2*x1+sqrt(5)*sin(t);x2=3+sqrt(2)*cos(t); y2=6*sin(t);plot(x1,y1,x2,y2); grid on; 作图发现4个解的大致位置,然后分别求解y1=fsolve('[(x(1)-2)^2+(x(2)-3+2*x(1))^2-5,2*(x(1)-3)^2+(x(2)/3)^2-4]',[1.5,2])y2=fsolve('[(x(1)-2)^2+(x(2)-3+2*x(1))^2-5,2*(x(1)-3)^2+(x(2)/3)^2-4]',[1.8,-2])y3=fsolve('[(x(1)-2)^2+(x(2)-3+2*x(1))^2-5,2*(x(1)-3)^2+(x(2)/3)^2-4]',[3.5,-5])y4=fsolve('[(x(1)-2)^2+(x(2)-3+2*x(1))^2-5,2*(x(1)-3)^2+(x(2)/3)^2-4]',[4,-4])Page84 Exercise 8(1)clear;fun=inline('x.^2.*sin(x.^2-x-2)');fplot(fun,[-2 2]);grid on; 作图观察x(1)=-2;x(3)=fminbnd(fun,-1,-0.5);x(5)=fminbnd(fun,1,2);fun2=inline('-x.^2.*sin(x.^2-x-2)');x(2)=fminbnd(fun2,-2,-1);x(4)=fminbnd(fun2,-0.5,0.5);x(6)=2feval(fun,x)答案: 以上x(1)(3)(5)是局部极小,x(2)(4)(6)是局部极大,从最后一句知道x(1)全局最小,x(2)最大。
《数学实验》课程教学大纲
3.线性方程组 的求解
4.随机数的生产和模拟
5.实验实例:循环比赛的名次和按年龄分组的种最优化问题实验
重点:学会一些常用函数的调用格式并学会自己动手编写函数
3. 《高等应用数学问题的MATLAB求解》.薛定宇,陈阳泉著.清华大学出版社,2004
4. 《MATLAB数学实验》.胡良剑,孙晓君编著.高等教育出版社,2006.6
执笔人:邓化宇
审核人:
院(系)负责人:
《数学实验》课程教学大纲
MathematicalExperiment
适用:本科四年制信息与计算科学专业(40学时左右)
一、课程的目的及任务
开设《数学实验》课的目的是在两周的时间里为学生介绍如何使用计算机的语言和方法去处理一些经典的数学问题,并提供一些实例以启发学生自己动手练习。进一步的提高要靠学生的兴趣和努力。
教学要点:
1.一元非线性方程数值求解
2.非线性方程组数值求解
3.方程符号求解
4.一元函数和多元函数无约束优化求解
5.线性规划
6.实验实例:购房贷款的利率和最短路问题
第五章 微分方程问题的计算机求解
重点:学会一些常用函数的调用格式并学会自己动手编写函数
教学要点:
1.常系数微分方程的计算机求解析解
2.微分方程问题的数值解法
二、课程的特点、要求及本课程与其它课程的联系
数学是科学技术人才科学素质的的重要组成部分,随着高科技与与计算技术的发展和普及,数学的重要性日益突出。“高技术本质上是一种数学技术”这一观点已越来越多地为人们所认同。学习计算机使用和开发是启迪学生创新意识和创新思维、锻炼创新能力、培养高层次人才的一条重要途径;也是激发学习欲望、培养主动探索、努力进取学风和团结协作精神的有力措施。
中国海洋大学本科生课程大纲
中国海洋大学本科生课程大纲课程名称数学实验Mathematics Experiments课程代码 0751********课程属性 工作技能 课时/学分64/2课程性质 必修 实践学时 64 责任教师 施心慧 课外学时 0 课程属性:公共基础/通识教育/学科基础/专业知识/工作技能,课程性质:必修、选修一、 课程介绍1.课程描述:数学实验是由于计算机技术和科学计算软件的迅猛发展应运而生的一门较新的数学课程,它改变了数学只靠纸和笔的传统形象,将实验的手段引入到数学的学习和研究中。
本课程为大学二年级数学院的学生开设。
它不是讲授新的数学知识,而是让学生利用已有的数学知识去解决一些经简化的实际问题。
大多数实验的一般过程是:对于给出的实际问题,建立数学模型、选择适当的数学方法、用科学计算软件MATLAB编程计算、对运算结果进行分析、给出结论。
本课程以MATLAB软件为主要的实验工具,采用以学生动手动脑为主,教师讲授和点评、小组讨论、报告为辅的教学方式。
通过本课程的学习,学生用数学解决实际问题的意识和能力可以得到强化和提高,更切实地体会到数学的用处,增加学习兴趣,提高创造力。
2.设计思路:本课程旨在训练用数学解决实际问题的能力。
实验内容的选取是基于学生具备MATLAB语言的初步编程能力、并学习了数学分析、高等代数、解析几何、运筹学基础(初步)、数学实验基础、常微分方程、数值分析或计算方法、概率论等数学课程的基础之上。
课程共分七个基础实验和一个综合实验依次进行。
七个基础实验是:MATLAB 基础知识复习、常微分方程(组)、数据建模——插值与拟合、古典密码学、图与网络- 1 -优化、动态规划、遗传算法。
基础实验涉及的数学内容较为单一、数学模型和求解方法较简单,是对“用数学”能力的基本训练。
综合实验以三人为一组进行,所涉及到的数学知识范围更广,建模和求解的难度更大。
综合实验的题目可以小组自拟或在任课教师拟定的题目中选择。
最全的Matlab经典教程免费下载
郭毓《MATLAB与SIMULINK基础》讲义--南京理工大学自动化系/thread-213-1-1.html
《Matlab6.0数学手册》pdf版(by Archan)/thread-202-1-2.html
《MATLAB实验仿真与算法设计》教学课件(谢红梅) /thread-201-1-2.html
冈萨雷斯《数字图像处理》中文版/thread-180-1-2.html
《精通MATLAB综合辅导与指南例程》讲义/thread-206-1-2.html
《科学计算与MATLAB语言》讲义 /thread-204-1-2.html
《MATLAB 程式设计与应用》讲义(by 张智星)/thread-203-1-2.html
sandy《Matlab与数值分析简明教程》/thread-320-1-1.html
满晓宇《战胜MATLAB必做练习50题》/thread-319-1-1.html
宋新山《Matlab在环境科学中的应用》/thread-303-1-1.html
苏金明《Matlab高级编程》 /thread-186-1-3.html
飞思科技《MATLAB6.5应用接口编程》 /thread-185-1-3.html
飞思科技《MATLAB6.5辅助小波分析与应用》 /thread-184-1-3.html
李国勇《智能控制及其Matlab实现》/thread-189-1-3.html
韩利竹《Matlab电子仿真与应用(第二版)》/thread-188-1-3.html
程卫国《Matlab5.3精要、编程及高级应用》/thread-187-1-3.html
20182019s课程教材信息
中国人民大学出版社 郭寿康、韩立余
课程 [040821外经贸会话] [040891纺织材料学] [040903非织造学] [042011服饰搭配, 047061服饰搭配与赏析] [045472纺织结构复合材料] [046121纺纱工艺学] [046682大提花织物设计] [046831微纳米纺织品与检测] [047101纺纱学(1)] [047431分析化学] [047501采购学] [047551针织面料与成形产品设计] [060036集成产品设计]
课程 [020761行政法与行政诉讼法学] []020881西方政治学说史 [020941西方行政学说史] [021001广告学] [022221整合营销传播] [022311知识产权法] [022672毛泽东思想和中国特色社会主义理论体 系概论Ⅱ] [022851传播学概论] [022881中外新闻传播史] [022901传播伦理与职业道德] [023051传播学研究方法] [023321纪录片理论与实务] [023511统计原理与应用] [023601新媒体文化产业] [024271微机原理] [024281认知与交互设计] [024321影视美学] [024351危机公关综合实务] [024381新媒体栏目策划与创意] [024441时尚传播史] [024461金融法(下)] [024571逻辑与批判性思维] [024711心灵与人工智能的哲学导论] [024741科学哲学导论] [024771法律文书] [024821国际投资法] [024841国际投资法基础] [040043针织学(1), 046131针织工艺学] [040091纺织科技史导论] [040232市场营销学] [040242国际贸易法]
中国人民大学出版社 胡百精 中国人民大学出版社 许鹏 中国纺织出版社 北京大学出版社 冯泽民 谷振诣 刘壮虎 中国人民大学出版社 朱大旗
MATLB数学实验(胡良剑)第十,十一十二章部分答案
第十章
第一题
f = [-3;4;-2;5];A =[1,1,3,-1;2,-3,1,-2]; b = [14;-2];Aeq = [4,-1,2,-1];beq = -2; i = 1:3;lb(i) = zeros(3,1); lb(4) = -inf; [x,feval] =linprog(f,A,b,Aeq,beq,lb) Optimization terminated.
rgeScale='off';end
if nargin<9,id=ones(size(f));end
if nargin<8,x=[];end
if nargin<7|isempty(ub),ub=inf*ones(size(f));end
if nargin<6|isempty(lb),lb=zeros(size(f));end
switching to medium-scale (line search).
> In fmincon at 260
Optimization terminated: first-order optimality measure less
than options.TolFun and maximum constraint violation is less
[tmpd,j]=min(W(i,V)); tmpj=V(j); for k=2:ndd
数学实验第四章函数和方程习题MATLAB编码(胡良剑版)
%p0,p1是初始值
%delta是给定允许误差
%m是p1-p0的误差估计?
%k是所需
%要的迭代次数?%y=f(p1)
% K=0,p0,p1,feval('f',p0),feval('f',p1)
% for k=1:max1
% if(f1==0)
% root=a;
% end
% if(f2==0)
% root=b;
% end
% if(f1*f2>0)
% disp('两端点函数值乘积大于0!');
% return;
% else
% tol=1;
% fun=diff(sym(f));
% fa=subs(sym(f),findsym(sym(f)),a);
% x=fminsearch(fun,[0 0]) %求极小值
% fun2=inline('-(x(2)^3/9+3*x(1)^2*x(2)+9*x(1)^2+x(2)^2+x(1)*x(2)+9)');
% x=fminsearch(fun2,[0 -5]) %求极大值
%第10题
% t=0:24;
% c=[15 14 14 14 14 15 16 18 20 22 23 25 28 ...
% x(2)=fminsearch(fun2,-2.5);
% x(4)=3;
% feval(fun,x)
%第8题(3)
% fun=inline('abs(x^3-x^2-x-2)');
% fplot(fun,[0 3]);grid on; %作图观察
基于Matlab的数值分析实验的设计r——以矩阵的奇异值分解为例
基于Matlab的数值分析实验的设计r——以矩阵的奇异值分解为例杨海涛【摘要】基于Matlab的数学实验,目的是为了提高普通高校学生的数学兴趣和动手能力,为培养解决实际数学问题的能力而设计的实验方法和步骤.利用矩阵的奇异值分解svd函数的数值算法设计了一种在实验课上用svd算法仿真的实验方案.【期刊名称】《内蒙古民族大学学报(自然科学版)》【年(卷),期】2016(031)006【总页数】4页(P465-468)【关键词】数学实验;数值计算;Matlab;奇异值分解【作者】杨海涛【作者单位】内蒙古民族大学数学学院,内蒙古通辽 028043【正文语种】中文【中图分类】O241.4利用计算机计算解决实际问题时需要建立数学模型,同时也需要相应数学软件的支持,为了更好的掌握数学软件这一有利的工具,本文通过Matlab中的矩阵的奇异值分解〔1〕svd函数设计了如何学习和掌握该函数的实验方法,从而可以进一步学习其它函数.目前,Matlab在计算仿真和图形图像的压缩处理〔2~4〕有着广泛的应用,希望这个实验方法能够给正在学习数学软件的学生起到抛砖引玉的作用.1.1 理论基础任意复长方矩阵的奇异值分解定理如下(Autonne-Eckart-Young定理):〔5〕设A为一个m×n的实矩阵(m≥n),则存在正交矩阵V和U,使得:其中∑=diag(σ1,σ2,…,σr),而且σ1≥σ2≥…≥σr≥0.奇异值分解在统计分析、信号与图像处理、系统理论和控制中有着广泛的应用. 1.2 SVD的数值算法〔6〕求解一般矩阵奇异值问题最常用的算法可分为两大类:一类是QR分解的方法;另一类是Jacobi旋转的方法.奇异值分解的QR分解法分为两个阶段:一是矩阵的二重对角化,即利用Householder变换把矩阵A变换成二重对角矩阵.这个过程就是追赶(chasing)过程〔7〕.二是QR分解,即保持二重对角矩阵形式不变,利用正交变换使对角线元素逐渐减小,使得矩阵接近对角矩阵.1.2.1 svd函数.对于matlab中的svd函数不是开源的,但通过调研得知该函数的算法就是采用propack中的lansvd函数的算法.其算法过程如下:对于任意的矩阵Am×n可以利用lanczos过程将其双对角化,即存在m阶正交矩阵Q和n阶正交矩阵W满足:其中为一双对角阵,且αi, βi∈R (i=1, 2, …, n),这里假设了m>n.上面的分解是逐步的利用Householder左右变换得到的.算法中利用了已知函数完成的双对角化过程.这里如果我们已经知道Bn的svd分解为其中Vn+1为n+1阶正交矩阵,的广义对角矩阵Un为n阶正交矩阵.则有令V=Qn+1Vn+1,其前r列为对应的左奇异向量,令U=WUn,其前r列为对应的右奇异向量.因此算法的关键在于计算矩阵Bn的svd分解.因为Bn是一个双对角矩阵,它的奇异值分解计算比较简单,运算量较小,利用Givens变换,算法复杂度仅为O(n).在我们这个程序中,直接调用了库函数svd 计算Bn的奇异值.1.2.2 svds函数.对于matlab中的svds函数,算法流程如下:(1)对于m×n阶矩阵A,构造(m+n)×(m+n)阶矩阵,其中O是稀疏矩阵. (2)利用eigs函数计算出B矩阵的特征值矩阵D,以及D对应的正交阵W. (3)根据输入参数,按要求输出最大特征值、最接近某参数的特征值或对特征值进行排序等.(4)根据设定的阀值,对矩阵进行筛选得到U矩阵,对矩阵筛选得到V矩阵,对D矩阵进行筛选得到S矩阵.(5)根据S中元素的大小对U、V、S进行排序.(6)输出结果.有如下J-W定理(Jordan-Wielandt定理):若σ1≥σ2≥…≥σp-1≥σp是矩阵Am×n的奇异值(其中,p=min{m,n}),则上述B矩阵具有2p个非零特征值-σ1,…,-σp,σp,…,σ1和|m-n|个零特征值,与非零特征值(±σj)相对应的特征向量为若m≠n,另有特征向量分别对应于m>n或m<n.svds函数的基本思想就是构造出这样的矩阵B,并且求出其特征值和特征向量,然后根据函数重构的条件(所需奇异值的数量、最大最小等),求出所需的奇异值和对应的向量.2.1 矩阵生成生成10个不同的m×n的矩阵,这里取m=120 , n=50,对这10个矩阵进行svd分解.矩阵生成的方法是先做出矩阵的奇异值,然后用QR分解做出10个120阶正交阵和10个50阶正交阵,再利用奇异值分解公式生成10个矩阵A.10个矩阵的奇异值(均为25个)如表1所示.2.2 仿真结果2.2.1 计算时间svd和svds两种方法对十个矩阵进行奇异值分解消耗的时间如表2,单位是秒.可见,svd函数运行时间一般要小于svds函数运行时间.2.2.2 最小奇异值计算准确率.通过对比最小特征值的计算准确度,将每种算法特征值的最小值与矩阵最开始生成时的最小特征值进行比较.比较结果两种方法对10个矩阵的svd分解运算最小特征值的误差均几乎为0,都在10-15次方的量级左右,见表3.其中delta_min_svd和delta_min_svds是matlab程序中两个变量值.这个误差已经和截断误差非常接近了,因此可以认为两种方法都能较好的进行奇异值计算.2.2.3 左右奇异值向量计算准确率.由于svd分解后U、V矩阵的不唯一性,这里我们只看生成的U、V的正交性是不是依然比较完好.通过计算每一列和它的转置相乘减去单位阵也都在10-15这样的量级,因此svd和svds对于正交性的保持都是非常好的.综上所述,通过上述设计的方法和步骤,完全可以掌握matlab中的M-文件和.mat文件的使用方法,并且能熟练掌握奇异值分解的理论知识和具体计算方法,才能为以后svd应用打下基础,这样的实验方法才有实际的意义,才是真正的综合设计实验.本文的宗旨在于引导学生如何用数学软件完成一个实验设计,同时掌握这个实验的目的、相关内容实验的方法和手段.〔1〕尹芳黎,杨雁莹,王传栋,等.矩阵奇异值分解及其在高维数据处理中的应用〔J〕.数学的实践与认识,2011,45(15): 171-176.〔2〕胡乡峰,卫金茂.基于奇异值分解(SVD)的图像压缩〔J〕.东北师大学报:自然科学版,2006,38(3):36-39.〔3〕曾超,张陈,徐永利.基于奇异值分解的图像压缩及其matlab实现〔J〕.科技信息,2010,14:484.〔4〕向培素.一种自适应AP算法的matlab实现〔J〕.西南民族大学学报:自然科学版,2014,40(16):877-882.〔5〕张贤达.矩阵分析与应用〔M〕.北京:清华大学出版社,2004.358-360. 〔6〕姜健飞,胡良剑,唐俭.数值分析及其matlab实验〔M〕.北京:科学出版社,2004.6,219-224.〔7〕山晓东,李园.一类不可约L-矩阵的预条件AOR迭代法〔J〕.湖北民族学院学报:自然科学版,2014,32(2):128-132.。
数学实验第二版答案
数学实验第二版答案【篇一:matlab数学实验习题全部答案(胡良剑)】ge20,ex1(5) 等于[exp(1),exp(2);exp(3),exp(4)](7) 3=1*3, 8=2*4(8) a为各列最小值,b为最小值所在的行号(10) 1=4,false, 2=3,false, 3=2, ture, 4=1,ture(11) 答案表明:编址第2元素满足不等式(30=20)和编址第4元素满足不等式(40=10)(12) 答案表明:编址第2行第1列元素满足不等式(30=20)和编址第2行第2列元素满足不等式(40=10)%page20, ex2(1)a, b, c的值尽管都是1,但数据类型分别为数值,字符,逻辑,注意a与c相等,但他们不等于b (2)double(fun)输出的分别是字符a,b,s,(,x,)的ascii码%page20,ex3 r=2;p=0.5;n=12;t=log(r)/n/log(1+0.01*p)t =11.5813%page20,ex4x=-2:0.05:2;f=x.^4-2.^x;[fmin,min_index]=min(f)fmin =-1.3907 %最小值min_index =54 %最小值点编址x(min_index)ans =0.6500 %最小值点[f1,x1_index]=min(abs(f)) %求近似根--绝对值最小的点f1 =0.0328x1_index =24x(x1_index)ans =x(x1_index)=[];f=x.^4-2.^x; %删去绝对值最小的点以求函数绝对值次小的点 [f2,x2_index]=min(abs(f)) %求另一近似根--函数绝对值次小的点f2 =0.0630x2_index =651x(x2_index)ans = 1.2500%page20,ex5z=magic(10)z =92 99 1 8 15 67 74 51 58 4098 80 7 14 16 73 55 57 64 414 81 88 20 22 54 56 63 70 4785 87 19 21 3 60 62 69 71 2886 93 25 2 9 61 68 75 52 3417 24 76 83 90 42 49 26 33 6523 5 82 89 91 48 30 32 39 6679 6 13 95 97 29 31 38 45 7210 12 94 96 78 35 37 44 46 5311 18 100 77 84 36 43 50 27 59sum(z)ans =505 505 505 505 505 505 505 505 505 505 sum(diag(z))ans =505z(:,2)/sqrt(3)ans =57.157746.188046.765450.229553.693613.85642.88683.46416.9282z(8,:)=z(8,:)+z(3,:)z =92 99 1 8 15 67 74 51 58 4098 80 7 14 16 73 55 57 64 414 81 88 20 22 54 56 63 70 4785 87 19 21 3 60 62 69 71 2886 93 25 2 9 61 68 75 52 3417 24 76 83 90 42 49 26 33 6523 5 82 89 91 48 30 32 39 66283 87 101 115 119 83 87 101 115 11910 12 94 96 78 35 37 44 46 5311 18 100 77 84 36 43 50 27 593%page 40 ex1先在编辑器窗口写下列m函数,保存为eg2_1.m function [xbar,s]=ex2_1(x)n=length(x);xbar=sum(x)/n;s=sqrt((sum(x.^2)-n*xbar^2)/(n-1));例如x=[81 70 65 51 76 66 90 87 61 77];[xbar,s]=ex2_1(x)xbar =72.4000s = 12.1124%page 40 ex2s=log(1);n=0;while s=100n=n+1;s=s+log(1+n);endm=n 计算结果m=37%page 40 ex3clear;f(1)=1;f(2)=1;k=2;x=0;e=1e-8; a=(1+sqrt(5))/2;while abs(x-a)ek=k+1;f(k)=f(k-1)+f(k-2); x=f(k)/f(k-1); end计算至k=21可满足精度%page 40 ex4clear;tic;s=0;for i=1:1000000s=s+sqrt(3)/2^i;ends,toctic;s=0;i=1;while i=1000000s=s+sqrt(3)/2^i;i=i+1;4ends,toctic;s=0;i=1:1000000;s=sqrt(3)*sum(1./2.^i); s,toc%page 40 ex5t=0:24;c=[15 14 14 14 14 15 16 18 20 22 23 25 28 ... 31 32 31 29 27 25 24 22 20 18 17 16]; plot(t,c)%page 40 ex6%(1)x=-2:0.1:2;y=x.^2.*sin(x.^2-x-2);plot(x,y) y=inline(x^2*sin(x^2-x-2));fplot(y,[-2 2]) %(2)参数方法t=linspace(0,2*pi,100);x=2*cos(t);y=3*sin(t); plot(x,y)%(3)x=-3:0.1:3;y=x;[x,y]=meshgrid(x,y);z=x.^2+y.^2;surf(x,y,z)%(4)x=-3:0.1:3;y=-3:0.1:13;[x,y]=meshgrid(x,y);z=x.^4+3*x.^2+y.^2-2*x-2*y-2*x.^2.*y+6; surf(x,y,z)%(5)t=0:0.01:2*pi;x=sin(t);y=cos(t);z=cos(2*t);plot3(x,y,z)theta=linspace(0,2*pi,50);fai=linspace(0,pi/2,20);[theta,fai]=meshgrid(theta,fai);x=2*sin(fai).*cos(theta);5【篇二:数学实验(matlab)课后习题答案】常见曲线的图形。
中国海洋大学本科生课程大纲-数学科学学院
四、教学进度
序号
第一章 1.1-1.3
第二章 2.1-2.2
专题 或主题
MATLAB 入门
MATLAB 编程与 作图
计划 课时
8
16
主要内容概述
MATLAB 桌面(安装与运行) 数据和变量 数组及其运算 关系与逻辑运算 建模试验:贷款利率模型 程序设计(分支、循环结构) MATLAB 作图(二维平面曲线图、三维空 间曲线图、三维空间曲面图) 动画、游戏制作(简介)
三、学习要求 要完成所有的课程任务,学生必须: (1)按时上课,认真听讲,积极参与课堂讨论、和上机练习。本课程将包含较多
的上机练习、小组讨论等课堂活动。 (2)认真对待并完成规定的上机任务。由于本课程实践性强,学生个人能力有差
别,因此上机任务极有可能无法在课堂上完成,因此需要学生能够利用课下时间,继 续完成布置的上机任务。这些任务能加深对课程内容的理解、促进独立查阅资料和解 决问题的能力培养。上机任务的完成情况,是平时成绩的关键组成部分。
中国海洋大学本科生课程大纲
课程名称 课程属性
数学实验基础 Foundation of Mathematics Experiments
工作技能
东华大学教材信息
022681 022691 022761 022891 022901 022931 022961 023051 023161
课程名称
应用写作 西方经济学
宪法学
教材名称
应用写作教程 经济学原理 微观经济分册
宪法学
民法学 国际公法
民法总论 国际法
刑事诉讼法学 律师与公证制度
合同法
刑事诉讼法 中国律师学
沈莘 杨正泉
韩松 本书编写组
本书编写组
罗豪才 湛中乐 陈岳
:张昕,李泉 谭云明 黄瑚 李培林
章培恒,骆玉明 艾尔 巴比 蒋承勇
第3页
面向专业
人文学院 法学 法学 法学 法学 法学 法学 法学
行政管理、公共关系等 法学
所选该课程的学生 行政管理、公共关系学
公管 行政 行政管理
行政 公共关系
全校
行政管理 行政11 公共关系
课程代码
023321
课程名称
纪录片理论与实务
教材名称
纪录片创作
出版社
中国人民大学出版社
作者
朱景和
面向专业
书库号
传播学、新闻学、广播电视 02-23
023461 023511 023601 023711 023741 023751 040063 040073 040091 040441 040651 040672 040681
多元统计分析
中国统计出版社
时间序列分析及应用:R语言
机械工业出版社
质量管理统计方法
中国统计出版社
中国精算师金融数学过关必做1000 中国石化出版社有限公
题:含历年真题
司
金融数学
中国财政经济出版社
东华大学MATLAB数学实验第二版答案(胡良剑)精编版
数学实验答案Chapter 1Page20,ex1(5) 等于[exp(1),exp(2);exp(3),exp(4)](7) 3=1*3, 8=2*4(8) a为各列最小值,b为最小值所在的行号(10) 1>=4,false, 2>=3,false, 3>=2, ture, 4>=1,ture(11) 答案表明:编址第2元素满足不等式(30>=20)和编址第4元素满足不等式(40>=10)(12) 答案表明:编址第2行第1列元素满足不等式(30>=20)和编址第2行第2列元素满足不等式(40>=10)Page20, ex2(1)a, b, c的值尽管都是1,但数据类型分别为数值,字符,逻辑,注意a与c相等,但他们不等于b(2)double(fun)输出的分别是字符a,b,s,(,x,)的ASCII码Page20,ex3>> r=2;p=0.5;n=12;>> T=log(r)/n/log(1+0.01*p)Page20,ex4>> x=-2:0.05:2;f=x.^4-2.^x;>> [fmin,min_index]=min(f)最小值最小值点编址>> x(min_index)ans =0.6500 最小值点>> [f1,x1_index]=min(abs(f)) 求近似根--绝对值最小的点f1 =0.0328x1_index =24>> x(x1_index)ans =-0.8500>> x(x1_index)=[];f=x.^4-2.^x; 删去绝对值最小的点以求函数绝对值次小的点>> [f2,x2_index]=min(abs(f)) 求另一近似根--函数绝对值次小的点f2 =0.0630x2_index =65>> x(x2_index)ans =1.2500>> z=magic(10)z =92 99 1 8 15 67 74 51 58 4098 80 7 14 16 73 55 57 64 414 81 88 20 22 54 56 63 70 4785 87 19 21 3 60 62 69 71 2886 93 25 2 9 61 68 75 52 3417 24 76 83 90 42 49 26 33 6523 5 82 89 91 48 30 32 39 6679 6 13 95 97 29 31 38 45 7210 12 94 96 78 35 37 44 46 5311 18 100 77 84 36 43 50 27 59>> sum(z)>> sum(diag(z))>> z(:,2)/sqrt(3)>> z(8,:)=z(8,:)+z(3,:)Chapter 2Page 45 ex1先在编辑器窗口写下列M函数,保存为eg2_1.m function [xbar,s]=ex2_1(x)n=length(x);xbar=sum(x)/n;s=sqrt((sum(x.^2)-n*xbar^2)/(n-1));例如>>x=[81 70 65 51 76 66 90 87 61 77];>>[xbar,s]=ex2_1(x)Page 45 ex2s=log(1);n=0;while s<=100n=n+1;s=s+log(1+n);endm=nPage 40 ex3clear;F(1)=1;F(2)=1;k=2;x=0;e=1e-8; a=(1+sqrt(5))/2;while abs(x-a)>ek=k+1;F(k)=F(k-1)+F(k-2); x=F(k)/F(k-1);enda,x,k计算至k=21可满足精度clear;tic;s=0;for i=1:1000000s=s+sqrt(3)/2^i;ends,toctic;s=0;i=1;while i<=1000000s=s+sqrt(3)/2^i;i=i+1;ends,toctic;s=0;i=1:1000000;s=sqrt(3)*sum(1./2.^i);s,tocPage 45 ex5t=0:24;c=[15 14 14 14 14 15 16 18 20 22 23 25 28 ...31 32 31 29 27 25 24 22 20 18 17 16];plot(t,c)Page 45 ex6(1)x=-2:0.1:2;y=x.^2.*sin(x.^2-x-2);plot(x,y)y=inline('x^2*sin(x^2-x-2)');fplot(y,[-2 2]) (2)参数方法t=linspace(0,2*pi,100);x=2*cos(t);y=3*sin(t); plot(x,y)(3)x=-3:0.1:3;y=x;[x,y]=meshgrid(x,y);z=x.^2+y.^2;surf(x,y,z)(4)x=-3:0.1:3;y=-3:0.1:13;[x,y]=meshgrid(x,y);z=x.^4+3*x.^2+y.^2-2*x-2*y-2*x.^2.*y+6;surf(x,y,z)(5)t=0:0.01:2*pi;x=sin(t);y=cos(t);z=cos(2*t);plot3(x,y,z)(6)theta=linspace(0,2*pi,50);fai=linspace(0,pi/2,20); [theta,fai]=meshgrid(theta,fai);x=2*sin(fai).*cos(theta);y=2*sin(fai).*sin(theta);z=2*cos(fai);surf(x,y,z)(7)x=linspace(0,pi,100);y1=sin(x);y2=sin(x).*sin(10*x);y3=-sin(x);plot(x,y1,x,y2,x,y3)page45, ex7x=-1.5:0.05:1.5;y=1.1*(x>1.1)+x.*(x<=1.1).*(x>=-1.1)-1.1*(x<-1.1);plot(x,y)page45,ex9clear;close;x=-2:0.1:2;y=x;[x,y]=meshgrid(x,y);a=0.5457;b=0.7575;p=a*exp(-0.75*y.^2-3.75*x.^2-1.5*x).*(x+y>1);p=p+b*exp(-y.^2-6*x.^2).*(x+y>-1).*(x+y<=1);p=p+a*exp(-0.75*y.^2-3.75*x.^2+1.5*x).*(x+y<=-1);mesh(x,y,p)page45, ex10lookfor lyapunovhelp lyap>> A=[1 2 3;4 5 6;7 8 0];C=[2 -5 -22;-5 -24 -56;-22 -56 -16];>> X=lyap(A,C)X =1.0000 -1.0000 -0.0000-1.0000 2.0000 1.0000-0.0000 1.0000 7.0000Chapter 3Page65 Ex1>> a=[1,2,3];b=[2,4,3];a./b,a.\b,a/b,a\bans =0.5000 0.5000 1.0000ans =2 2 1ans =0.6552 一元方程组x[2,4,3]=[1,2,3]的近似解ans =0 0 00 0 00.6667 1.3333 1.0000矩阵方程[1,2,3][x11,x12,x13;x21,x22,x23;x31,x32,x33]=[2,4,3]的特解Page65 Ex 2(1)>> A=[4 1 -1;3 2 -6;1 -5 3];b=[9;-2;1];>> rank(A), rank([A,b]) [A,b]为增广矩阵ans =3ans =3 可见方程组唯一解>> x=A\bx =2.38301.48942.0213(2)>> A=[4 -3 3;3 2 -6;1 -5 3];b=[-1;-2;1];>> rank(A), rank([A,b])ans =3ans =3 可见方程组唯一解>> x=A\bx =-0.4706-0.2941(3)>> A=[4 1;3 2;1 -5];b=[1;1;1];>> rank(A), rank([A,b])ans =2ans =3 可见方程组无解>> x=A\bx =0.3311-0.1219 最小二乘近似解(4)>> a=[2,1,-1,1;1,2,1,-1;1,1,2,1];b=[1 2 3]';%注意b的写法>> rank(a),rank([a,b])ans =3ans =3 rank(a)==rank([a,b])<4说明有无穷多解>> a\bans =110 一个特解Page65 Ex3>> a=[2,1,-1,1;1,2,1,-1;1,1,2,1];b=[1,2,3]';>> x=null(a),x0=a\bx =-0.62550.6255-0.20850.4170x0 =11通解kx+x0Page65 Ex 4>> x0=[0.2 0.8]';a=[0.99 0.05;0.01 0.95];>> x1=a*x, x2=a^2*x, x10=a^10*x>> x=x0;for i=1:1000,x=a*x;end,xx =0.83330.1667>> x0=[0.8 0.2]';>> x=x0;for i=1:1000,x=a*x;end,xx =0.83330.1667>> [v,e]=eig(a)v =0.9806 -0.70710.1961 0.7071e =1.0000 00 0.9400>> v(:,1)./xans =1.17671.1767 成比例,说明x是最大特征值对应的特征向量Page65 Ex5用到公式(3.11)(3.12)>> B=[6,2,1;2.25,1,0.2;3,0.2,1.8];x=[25 5 20]'; >> C=B/diag(x)C =0.2400 0.4000 0.05000.0900 0.2000 0.01000.1200 0.0400 0.0900>> A=eye(3,3)-CA =0.7600 -0.4000 -0.0500-0.0900 0.8000 -0.0100-0.1200 -0.0400 0.9100>> D=[17 17 17]';x=A\Dx =37.569625.786224.7690Page65 Ex 6(1)>> a=[4 1 -1;3 2 -6;1 -5 3];det(a),inv(a),[v,d]=eig(a) ans =-94ans =0.2553 -0.0213 0.04260.1596 -0.1383 -0.22340.1809 -0.2234 -0.0532v =0.0185 -0.9009 -0.3066-0.7693 -0.1240 -0.7248-0.6386 -0.4158 0.6170d =-3.0527 0 00 3.6760 00 0 8.3766(2)>> a=[1 1 -1;0 2 -1;-1 2 0];det(a),inv(a),[v,d]=eig(a) ans =1ans =2.0000 -2.0000 1.00001.0000 -1.0000 1.00002.0000 -3.0000 2.0000v =-0.5773 0.5774 + 0.0000i 0.5774 - 0.0000i-0.5773 0.5774 0.5774-0.5774 0.5773 - 0.0000i 0.5773 + 0.0000id =1.0000 0 00 1.0000 + 0.0000i 00 0 1.0000 - 0.0000i(3)>> A=[5 7 6 5;7 10 8 7;6 8 10 9;5 7 9 10]A =5 76 57 10 8 76 8 10 95 7 9 10>> det(A),inv(A), [v,d]=eig(A)ans =1ans =68.0000 -41.0000 -17.0000 10.0000-41.0000 25.0000 10.0000 -6.0000-17.0000 10.0000 5.0000 -3.000010.0000 -6.0000 -3.0000 2.0000v =0.8304 0.0933 0.3963 0.3803-0.5016 -0.3017 0.6149 0.5286-0.2086 0.7603 -0.2716 0.55200.1237 -0.5676 -0.6254 0.5209d =0.0102 0 0 00 0.8431 0 00 0 3.8581 00 0 0 30.2887(4)(以n=5为例)方法一(三个for)n=5;for i=1:n, a(i,i)=5;endfor i=1:(n-1),a(i,i+1)=6;endfor i=1:(n-1),a(i+1,i)=1;enda方法二(一个for)n=5;a=zeros(n,n);a(1,1:2)=[5 6];for i=2:(n-1),a(i,[i-1,i,i+1])=[1 5 6];enda(n,[n-1 n])=[1 5];a方法三(不用for)n=5;a=diag(5*ones(n,1));b=diag(6*ones(n-1,1));c=diag(ones(n-1,1));a=a+[zeros(n-1,1),b;zeros(1,n)]+[zeros(1,n);c,zeros(n-1,1)] 下列计算>> det(a)ans =665>> inv(a)ans =0.3173 -0.5865 1.0286 -1.6241 1.9489-0.0977 0.4887 -0.8571 1.3534 -1.62410.0286 -0.1429 0.5429 -0.8571 1.0286-0.0075 0.0376 -0.1429 0.4887 -0.58650.0015 -0.0075 0.0286 -0.0977 0.3173>> [v,d]=eig(a)v =-0.7843 -0.7843 -0.9237 0.9860 -0.92370.5546 -0.5546 -0.3771 -0.0000 0.3771-0.2614 -0.2614 0.0000 -0.1643 0.00000.0924 -0.0924 0.0628 -0.0000 -0.0628-0.0218 -0.0218 0.0257 0.0274 0.0257d =0.7574 0 0 0 00 9.2426 0 0 00 0 7.4495 0 00 0 0 5.0000 00 0 0 0 2.5505Page65 Ex 7(1)>> a=[4 1 -1;3 2 -6;1 -5 3];[v,d]=eig(a)v =0.0185 -0.9009 -0.3066-0.7693 -0.1240 -0.7248-0.6386 -0.4158 0.6170d =-3.0527 0 00 3.6760 00 0 8.3766>> det(v)ans =-0.9255 %v行列式正常, 特征向量线性相关,可对角化>> inv(v)*a*v 验算ans =-3.0527 0.0000 -0.00000.0000 3.6760 -0.0000-0.0000 -0.0000 8.3766>> [v2,d2]=jordan(a) 也可用jordanv2 =0.0798 0.0076 0.91270.1886 -0.3141 0.1256-0.1605 -0.2607 0.4213 特征向量不同d2 =8.3766 0 00 -3.0527 - 0.0000i 00 0 3.6760 + 0.0000i>> v2\a*v2ans =8.3766 0 0.00000.0000 -3.0527 0.00000.0000 0.0000 3.6760>> v(:,1)./v2(:,2) 对应相同特征值的特征向量成比例ans =2.44912.44912.4491(2)>> a=[1 1 -1;0 2 -1;-1 2 0];[v,d]=eig(a)v =-0.5773 0.5774 + 0.0000i 0.5774 - 0.0000i-0.5773 0.5774 0.5774-0.5774 0.5773 - 0.0000i 0.5773 + 0.0000id =1.0000 0 00 1.0000 + 0.0000i 00 0 1.0000 - 0.0000i>> det(v)ans =-5.0566e-028 -5.1918e-017i v的行列式接近0, 特征向量线性相关,不可对角化>> [v,d]=jordan(a)v =1 0 11 0 01 -1 0d =1 1 00 1 10 0 1 jordan标准形不是对角的,所以不可对角化(3)>> A=[5 7 6 5;7 10 8 7;6 8 10 9;5 7 9 10]A =5 76 57 10 8 76 8 10 95 7 9 10>> [v,d]=eig(A)v =0.8304 0.0933 0.3963 0.3803-0.5016 -0.3017 0.6149 0.5286-0.2086 0.7603 -0.2716 0.55200.1237 -0.5676 -0.6254 0.5209d =0.0102 0 0 00 0.8431 0 00 0 3.8581 00 0 0 30.2887>> inv(v)*A*vans =0.0102 0.0000 -0.0000 0.00000.0000 0.8431 -0.0000 -0.0000-0.0000 0.0000 3.8581 -0.0000-0.0000 -0.0000 0 30.2887本题用jordan不行, 原因未知(4)参考6(4)和7(1)Page65 Exercise 8只有(3)对称, 且特征值全部大于零, 所以是正定矩阵. Page65 Exercise 9(1)>> a=[4 -3 1 3;2 -1 3 5;1 -1 -1 -1;3 -2 3 4;7 -6 -7 0]>> rank(a)ans =3>> rank(a(1:3,:))ans =2>> rank(a([1 2 4],:)) 1,2,4行为最大无关组ans =3>> b=a([1 2 4],:)';c=a([3 5],:)';>> b\c 线性表示的系数ans =0.5000 5.0000-0.5000 1.00000 -5.0000Page65 Exercise 10>> a=[1 -2 2;-2 -2 4;2 4 -2]>> [v,d]=eig(a)v =0.3333 0.9339 -0.12930.6667 -0.3304 -0.6681-0.6667 0.1365 -0.7327d =-7.0000 0 00 2.0000 00 0 2.0000>> v'*vans =1.0000 0.0000 0.00000.0000 1.0000 00.0000 0 1.0000 v确实是正交矩阵Page65 Exercise 11设经过6个电阻的电流分别为i1, ..., i6. 列方程组如下20-2i1=a; 5-3i2=c; a-3i3=c; a-4i4=b; c-5i5=b; b-3i6=0;i1=i3+i4;i5=i2+i3;i6=i4+i5;计算如下>> A=[1 0 0 2 0 0 0 0 0;0 0 1 0 3 0 0 0 0;1 0 -1 0 0 -3 0 0 0; 1 -1 0 0 0 0 -4 0 0;0 -1 1 0 0 0 0 -5 0;0 1 0 0 0 0 0 0 -3; 0 0 0 1 0 -1 -1 0 0;0 0 0 0 -1 -1 0 1 0;0 0 0 0 0 0 -1 -1 1];>>b=[20 5 0 0 0 0 0 0 0]'; A\bans =13.34536.44018.54203.3274-1.18071.60111.72630.42042.1467Page65 Exercise 12>> A=[1 2 3;4 5 6;7 8 0];>> left=sum(eig(A)), right=sum(trace(A))left =6.0000right =6>> left=prod(eig(A)), right=det(A) 原题有错, (-1)^n应删去left =27.0000right =27>> fA=(A-p(1)*eye(3,3))*(A-p(2)*eye(3,3))*(A-p(3)*eye(3,3)) fA =1.0e-012 *0.0853 0.1421 0.02840.1421 0.1421 0-0.0568 -0.1137 0.1705>> norm(fA) f(A)范数接近0ans =2.9536e-013Chapter 4Page84 Exercise 1(1)roots([1 1 1])(2)roots([3 0 -4 0 2 -1])(3)p=zeros(1,24);p([1 17 18 22])=[5 -6 8 -5];roots(p)(4)p1=[2 3];p2=conv(p1, p1);p3=conv(p1, p2);p3(end)=p3(end)-4; %原p3最后一个分量-4roots(p3)Page84 Exercise 2fun=inline('x*log(sqrt(x^2-1)+x)-sqrt(x^2-1)-0.5*x');fzero(fun,2)Page84 Exercise 3fun=inline('x^4-2^x');fplot(fun,[-2 2]);grid on;fzero(fun,-1),fzero(fun,1),fminbnd(fun,0.5,1.5)Page84 Exercise 4fun=inline('x*sin(1/x)','x');fplot(fun, [-0.1 0.1]);x=zeros(1,10);for i=1:10, x(i)=fzero(fun,(i-0.5)*0.01);end;x=[x,-x]Page84 Exercise 5fun=inline('[9*x(1)^2+36*x(2)^2+4*x(3)^2-36;x(1)^2-2*x(2)^2-20*x(3);16*x(1)-x(1)^3-2*x(2)^ 2-16*x(3)^2]','x');[a,b,c]=fsolve(fun,[0 0 0])Page84 Exercise 6fun=@(x)[x(1)-0.7*sin(x(1))-0.2*cos(x(2)),x(2)-0.7*cos(x(1))+0.2*sin(x(2))];[a,b,c]=fsolve(fun,[0.5 0.5])Page84 Exercise 7clear; close; t=0:pi/100:2*pi;x1=2+sqrt(5)*cos(t); y1=3-2*x1+sqrt(5)*sin(t);x2=3+sqrt(2)*cos(t); y2=6*sin(t);plot(x1,y1,x2,y2); grid on; 作图发现4个解的大致位置,然后分别求解y1=fsolve('[(x(1)-2)^2+(x(2)-3+2*x(1))^2-5,2*(x(1)-3)^2+(x(2)/3)^2-4]',[1.5,2])y2=fsolve('[(x(1)-2)^2+(x(2)-3+2*x(1))^2-5,2*(x(1)-3)^2+(x(2)/3)^2-4]',[1.8,-2])y3=fsolve('[(x(1)-2)^2+(x(2)-3+2*x(1))^2-5,2*(x(1)-3)^2+(x(2)/3)^2-4]',[3.5,-5])y4=fsolve('[(x(1)-2)^2+(x(2)-3+2*x(1))^2-5,2*(x(1)-3)^2+(x(2)/3)^2-4]',[4,-4])Page84 Exercise 8(1)clear;fun=inline('x.^2.*sin(x.^2-x-2)');fplot(fun,[-2 2]);grid on; 作图观察x(1)=-2;x(3)=fminbnd(fun,-1,-0.5);x(5)=fminbnd(fun,1,2);fun2=inline('-x.^2.*sin(x.^2-x-2)');x(2)=fminbnd(fun2,-2,-1);x(4)=fminbnd(fun2,-0.5,0.5);x(6)=2feval(fun,x)答案: 以上x(1)(3)(5)是局部极小,x(2)(4)(6)是局部极大,从最后一句知道x(1)全局最小,x(2)最大。
计算方法1-2
e ( x*)
( er ( x ) r )
为
x的相对误差限。
3.有效数字 m 设 x 10 (0.a1 a 2 a n a p ) (a1 0, p ) 1 mn 若 (1.1) x * x 10 2
则说
x 具有n位有效数字,分别是
若
n p ,则称
10 0.04 10 0.1 10 0.00000000004 10
8 9 9
8
3.尽量避免相近数相减 例如,当x很大时,应
,
x1
x
1 x x1
1 1 1 x x 1 x ( x 1) 当x接近于0时,应
1 cos x sin x x 或tg sin x 1 cos x 2
f(x)=0根或f(x)零点,当f(x)复杂时,很难求 (找近似有效简单方法)。
§2.1 区间二分法
理 论 : f(x) ∈ C[a,b],单调, f(a)f(b)<0
f(x)=0在(a,b)有惟一根。
方程求根的步骤:(1)求根的隔离区间;(2)将根精确化 根的隔离区间求法:画草图;多项式函数交点横坐标; 试算。
程序运行 结果:
基本要求:
1.熟悉计算方法在解决实际问题中所处的地位,熟 悉计算方法是以计算机为工具求近似解的数值方 法; 2.熟悉绝对误差(限),相对误差(限)及有效数 字概念; 3.熟悉公式(1.2)--(1.9); 4.熟悉选用算法应遵循的原则;
作业:
第一章 8, 10
第二章 方程求根
§1.2 误差基础知识
一 .误差来源(分类)
1. 模型误差。 2. 观测误差。 3. 截断误差,如
x x sin x x ...... , 3! 5!
《数学实验》教学大纲
《数学实验》教学大纲课程编号:10107014学时:18学分:1课程类别:限制性选修课面向对象:数学与应用数学专业本科学生课程英文名称:Experiments in Mathematics一、课程的任务和目的任务:通过实验,使学生熟练掌握常用的数学实验方法,培养和提高应用计算机和相关软件进行科学与工程计算的能力,为以后的学习及应用打下良好的基础。
目的:使学生掌握数学实验的基本思想和方法,即不把数学看成先验的逻辑体系,而是把它视为一门“实验科学”,从问题出发,通过学习Matlab或Mathematica数学软件,借助计算机,学生亲自设计和动手,体验解决问题的过程,从实验中去学习、探索和发现数学的规律。
二、课程教学内容与要求(一)MATLAB桌面操作、矩阵与数组函数1.教学内容MATLAB桌面、矩阵与数组的建立及运算、Matlab常用数学函数、字符串、输入、输出。
2.基本要求掌握Matlab语言概况及基本知识。
3.重点和难点重点是语言基本知识的掌握,难点是Matlab语言与其他语言的差异。
(二) MATLAB程序设计与绘图1.教学内容MATLAB表达式、选择结构、循环结构、M文件(函数文件的定义)、全局变量和局部变量、数据和函数的可视化。
2.基本要求掌握MATLAB语言程序的设计与编写,能编程解决一些数学问题及绘图。
3.重点和难点重点是程序的设计与编写、错误的诊断和修正,难点是程序的设计、错误的诊断和修正。
(三)插值与拟合1.教学内容拉格朗日插值、分段线性插值、三次样条插值、曲线拟合的线性最小二乘法。
2.基本要求掌握用MATLAB计算拉格朗日、分段线性、三次样条插值的方法,能编写用线性最小乘法来作曲线拟合的程序。
3.重点和难点重点是一维插值和用多项式来拟合函数,难点是插值与拟合的源程序编写。
(四)应用微积分1.教学内容数值微积分MATLAB命令、计算实验:数值微积分、建模实验:奶油蛋糕。
2.基本要求掌握用MATLAB进行数值计算的方法,能编程计算数值导数和数值积分。
MATLAB习题参考答案(胡良剑
第一章 MATLAB 入门4、求近似解解:>> x=-2:0.05:2;y=x.^4-2.^x两个近似解:y1=f(-0.85)= -0.0328; y2=f(1.250)= 0.0630第二章 MATLAB 编程与作图1、 设x 是数组,求均值和方差解:函数文件如下:function [xx,s]=func1(x)n=length(x);xx=sum(x)/n;s=sqrt((sum(x.^2)-n*xx^2)/(n-1));命令窗口:>> x=[1 2 3 4 5];[xx,s]=func1(x)2、求满足的最小m 值 100)1ln(0>+∑=m n n s=0;n=0;while(s<=100)s=s+log(1+n);n=n+1;endn,s3、用循环语句形成Fibonacci 数列,....4,3,,12121=+===−−k F F F F F k k k 。
并验证极限2511+→−k k F F (提示:计算至两边误差小于精度1e-8为止) 解: 求Fibonacci 数列的函数文件:function f=fun(n)if n<=2f=1;elsef=fun(n-1)+fun(n-2);end验证极限的函数文件:function [k,a]=funTest(e)a=abs(1-(1+sqrt(5))/2);k=2;while(a>e)k=k+1;a=abs(fun(k)/fun(k-1)-(1+sqrt(5))/2);end命令行:>> [k,a]=funTest(10^-8)k =21a =9.7719e-009或者M 文件如下:clear; F(1)=1;F(2)=1;k=2;x=0;e=1e-8; a=(1+sqrt(5))/2;while abs(x-a)>ek=k+1; F(k)=F(k-1)+F(k-2); x=F(k)/F(k-1);enda,x,k4、分别用for 和while 循环结构编写程序,求出∑==610123i i K ,并考虑一种避免循环语句的程序设计,比较各种算法的运行时间。
matlab数学实验胡良剑第五章
第五章%Exercise 1x=[0 4 10 12 15 22 28 34 40];y=[0 1 3 6 8 9 5 3 0];trapz(x,y)%Exercise 2x=[0 4 10 12 15 22 28 34 40];y=[0 1 3 6 8 9 5 3 0];diff(y)./diff(x)%Exercise 3xa=-1:0.1:1;ya=0:0.1:2;[x,y]=meshgrid(xa,ya);z=x.*exp(-x.^2 -y.^3);[px,py] = gradient(z,xa,ya);px%Exercise 4t=0:0.01:1.5;x=log(cos(t));y=cos(t)-t.*sin(t);dydx=gradient(y,x)plot(x,dydx) %dydx函数图,作图观察x=-1时,dydx的值约0.9 %以下是更精确的编程计算方法[x_1,id]=min(abs(x-(-1)));%找最接近x=-1的点,id为这个点的下标dydx(id)%Exercise 5(1)fun=@(x)1/sqrt(2*pi)*exp(-x.^2/2);quadl(fun,0,1)或用trapzx=linspace(0,1,100);y=1/sqrt(2*pi)*exp(-x.^2/2);trapz(x,y)%Exercise 5(2)fun=inline('exp(2*x).*cos(x).^3');quadl(fun,0,2*pi)或用trapzx=linspace(0,2*pi,100);y=exp(2*x).*cos(x).^3;trapz(x,y)%Exercise 5(3)fun=@(x)x.*log(x.^4).*asin(1./x.^2);quadl(fun,1,3)或用trapzx=1:0.01:3;y=feval(fun,x);trapz(x,y)%Exercise 5(4)fun=@(x)sin(x)./x;quadl(fun,1e-10,1) %注意由于下限为0,被积函数没有意义,用很小的1e-10代替%Exercise 5(5)fun=inline('x.^(-x)','x');quadl(fun,1e-10,1) %注意由于下限为0,被积函数没有意义,用很小的1e-10代替%Exercise 5(6)fun=inline('sqrt(1+r.^2.*sin(th))','r','th');dblquad(fun,0,1,0,2*pi)%Exercise 5(7)%先在Editer窗口建立90页函数dblquad2,再在Command窗口clear;fun=@(x,y)1+x+y.^2;clo=@(x)-sqrt(2*x-x.^2);dhi=@(x)sqrt(2*x-x.^2);dblquad2(fun,0,2,clo,dhi,100)%Exercise 6t=linspace(0,2*pi,100);x=2*cos(t);y=3*sin(t);dx=gradient(x,t);dy=gradient(y,t);f=sqrt(dx.^2+dy.^2);trapz(t,f)%Exercise 6另一解法%先写参数方程x=2*cos(t);y=3*sin(t);%计算x'(t)=-2*sin(t),y'(t)=3*cos(t)%4*sin(t)^2+9*cos(t)^2=4+5*cos(t)^2fun=@(t)sqrt(4+5*cos(t).^2);quadl(fun,0,2*pi)%Exercise 7%先算出z的梯度dz/dx=(1-2*x^2)*exp(-x^2-y^2),dz/dy=(1-2*y^2)*exp(-x^2-y^2);%根据曲面面积公式fun=@(x,y)sqrt(1+((1-2*x.^2).*exp(-x.^2-y.^2)).^2+((1-2*y.^2)*exp(-x.^2-y.^2)).^2);dblquad(fun,-1,1,0,2)%或者用下列纯粹离散化解法xa=linspace(-1,1);ya=linspace(0,2);[x,y]=meshgrid(xa,ya);z=x.*exp(-x.^2-y.^2);[zx,zy]=gradient(z,xa,ya);f=sqrt(1+zx.^2+zy.^2);s=0;for i=2:length(xa)for j=2:length(ya)s=s+(xa(i)-xa(i-1))*(ya(j)-ya(j-1))*(f(i,j)+f(i-1,j)+f(i,j-1)+f(i-1,j-1))/4;%每个近似长方体高用四顶点平均值endends%Exercise 8funl=inline('-(-x).^0.2.*cos(x)');funr=inline('x.^0.2.*cos(x)');quadl(funl,-1,0)+quadl(funr,0,1)%Exercise 9 (以I32为例)fun=@(x)abs(sin(x));h=0.1;x=0:h:32*pi;y=feval(fun,x);t1=trapz(x,y)h=pi;x=0:h:32*pi;y=feval(fun,x);t2=trapz(x,y)%步长与周期一致,结果失真q1=quad(fun,0,32*pi)q2=quadl(fun,0,32*pi)%Exercise 10(1)先在Editer窗口建立88页函数deriv,再在Command窗口fun=inline('x.^2.*sin(x.^2+3*x-4)','x');deriv(fun,[1.3 1.5],0.1,1e-3) %取0.1为初始步长%注:书后习题答案错,1.3处导数应为2.4177,1.5处导数应为-11.3330%Exercise 10(2)%先在程序编辑器,写下列函数,保存为ex5_10_2ffunction d=ex5_10_2f(fname,a,h0,e)h=h0;d=(fname(a+h)-2*fname(a)+fname(a-h))/(h*h);d0=d+2*e;while abs(d-d0)>ed0=d;h0=h;h=h0/2;d=(fname(a+h)-2*fname(a)+fname(a-h))/(h*h);end%再在指令窗口执行fun=@(x)x.^2*sin(x.^2-x-2);d=ex5_10_2f(fun,1.4,0.1,1e-3)%Exercise 11%提示:f上升时,f'>0;f下降时,f'<0; f极值,f'=0.%Exercise 12在程序编辑器,写下列函数,保存为ex5_12function I=ex5_12(fname,a,b,n)h=(b-a)/n;x=a:h:b;f=fname(x);I=f(1)+f(n+1);for i=2:nif i-2*floor(i/2)==0I=I+4*f(i);elseI=I+2*f(i);endendI=h/3*I;%再在指令窗口执行ex5_12(inline('1/sqrt(2*pi)*exp(-x.^2/2)'),0,1,50) %注:原题n=5改为偶数n=50%更加符合Matlab风格的编程ex5_12function I=ex5_12f(fname,a,b,n)x=linspace(a,b,n+1);f=fname(x);I=(b-a)/n/3*(f(1)+f(n+1)+2*sum(f(3:2:n))+4*sum(f(2:2:n)));%Exercise 13fun=inline('5400*v./(8.276*v.^2+2000)','v');quadl(fun,15,30)%Exercise 14重心不超过凳边沿。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
数学实验答案Chapter 1Page20,ex1(5)等于[exp(1),exp(2);exp(3),exp(4)](7) 3=1*3, 8=2*4(8) a为各列最小值,b为最小值所在的行号(10) 1>=4,false, 2>=3,false, 3>=2, ture, 4>=1,ture(11)答案表明:编址第2元素满足不等式(30>=20)和编址第4元素满足不等式(40>=10)(12)答案表明:编址第2行第1列元素满足不等式(30>=20)和编址第2行第2列元素满足不等式(40>=10)Page20, ex2(1)a, b, c的值尽管都是1,但数据类型分别为数值,字符,逻辑,注意a与c相等,但他们不等于b(2)double(fun)输出的分别是字符a,b,s,(,x,)的ASCII码Page20,ex3>> r=2;p=0.5;n=12;>> T=log(r)/n/log(1+0.01*p)Page20,ex4>> x=-2:0.05:2;f=x.^4-2.^x;>> [fmin,min_index]=min(f)最小值最小值点编址>> x(min_index)ans =0.6500最小值点>> [f1,x1_index]=min(abs(f))求近似根--绝对值最小的点f1 =0.0328x1_index =24>> x(x1_index)ans =-0.8500>> x(x1_index)=[];f=x.^4-2.^x;删去绝对值最小的点以求函数绝对值次小的点>> [f2,x2_index]=min(abs(f))求另一近似根--函数绝对值次小的点f2 =0.0630x2_index =65>> x(x2_index)ans =1.2500Page20,ex5>> z=magic(10)z =92 99 1 8 15 67 74 51 58 4098 80 7 14 16 73 55 57 64 414 81 88 20 22 54 56 63 70 4785 87 19 21 3 60 62 69 71 2886 93 25 2 9 61 68 75 52 3417 24 76 83 90 42 49 26 33 6523 5 82 89 91 48 30 32 39 6679 6 13 95 97 29 31 38 45 7210 12 94 96 78 35 37 44 46 5311 18 100 77 84 36 43 50 27 59>> sum(z)>> sum(diag(z))>> z(:,2)/sqrt(3)>> z(8,:)=z(8,:)+z(3,:)Chapter 2Page 45 ex1先在编辑器窗口写下列M函数,保存为eg2_1.m function [xbar,s]=ex2_1(x)n=length(x);xbar=sum(x)/n;s=sqrt((sum(x.^2)-n*xbar^2)/(n-1));例如>>x=[81 70 65 51 76 66 90 87 61 77]; >>[xbar,s]=ex2_1(x)Page 45 ex2s=log(1);n=0;while s<=100n=n+1;s=s+log(1+n);endm=nPage 40 ex3clear;F(1)=1;F(2)=1;k=2;x=0;e=1e-8; a=(1+sqrt(5))/2;while abs(x-a)>ek=k+1;F(k)=F(k-1)+F(k-2); x=F(k)/F(k-1); enda,x,k计算至k=21可满足精度Page 45 ex4clear;tic;s=0;for i=1:1000000s=s+sqrt(3)/2^i;ends,toctic;s=0;i=1;while i<=1000000s=s+sqrt(3)/2^i;i=i+1;ends,toctic;s=0;i=1:1000000;s=sqrt(3)*sum(1./2.^i);s,tocPage 45 ex5t=0:24;c=[15 14 14 14 14 15 16 18 20 22 23 25 28 ...31 32 31 29 27 25 24 22 20 18 17 16];plot(t,c)Page 45 ex6(1)x=-2:0.1:2;y=x.^2.*sin(x.^2-x-2);plot(x,y)y=inline('x^2*sin(x^2-x-2)');fplot(y,[-2 2])(2)参数方法t=linspace(0,2*pi,100);x=2*cos(t);y=3*sin(t); plot(x,y)(3)x=-3:0.1:3;y=x;[x,y]=meshgrid(x,y);z=x.^2+y.^2;surf(x,y,z)(4)x=-3:0.1:3;y=-3:0.1:13;[x,y]=meshgrid(x,y);z=x.^4+3*x.^2+y.^2-2*x-2*y-2*x.^2.*y+6;surf(x,y,z)(5)t=0:0.01:2*pi;x=sin(t);y=cos(t);z=cos(2*t);plot3(x,y,z)(6)theta=linspace(0,2*pi,50);fai=linspace(0,pi/2,20);[theta,fai]=meshgrid(theta,fai);x=2*sin(fai).*cos(theta);y=2*sin(fai).*sin(theta);z=2*cos(fai);surf(x,y,z)(7)x=linspace(0,pi,100);y1=sin(x);y2=sin(x).*sin(10*x);y3=-sin(x);plot(x,y1,x,y2,x,y3)page45, ex7x=-1.5:0.05:1.5;y=1.1*(x>1.1)+x.*(x<=1.1).*(x>=-1.1)-1.1*(x<-1.1); plot(x,y)page45,ex9clear;close;x=-2:0.1:2;y=x;[x,y]=meshgrid(x,y);a=0.5457;b=0.7575;p=a*exp(-0.75*y.^2-3.75*x.^2-1.5*x).*(x+y>1);p=p+b*exp(-y.^2-6*x.^2).*(x+y>-1).*(x+y<=1);p=p+a*exp(-0.75*y.^2-3.75*x.^2+1.5*x).*(x+y<=-1); mesh(x,y,p)page45, ex10lookfor lyapunovhelp lyap>> A=[1 2 3;4 5 6;7 8 0];C=[2 -5 -22;-5 -24 -56;-22 -56 -16];>> X=lyap(A,C) X =1.0000 -1.0000 -0.0000-1.0000 2.0000 1.0000-0.0000 1.0000 7.0000Chapter 3Page65 Ex1>> a=[1,2,3];b=[2,4,3];a./b,a.\b,a/b,a\bans =0.5000 0.5000 1.0000ans =2 2 1ans =0.6552一元方程组x[2,4,3]=[1,2,3]的近似解ans =0 0 00 0 00.6667 1.3333 1.0000矩阵方程[1,2,3][x11,x12,x13;x21,x22,x23;x31,x32,x33]=[2,4,3]的特解Page65 Ex 2(1)>> A=[4 1 -1;3 2 -6;1 -5 3];b=[9;-2;1];>> rank(A), rank([A,b]) [A,b]为增广矩阵ans =3ans =3可见方程组唯一解>> x=A\bx =2.38301.48942.0213(2)>> A=[4 -3 3;3 2 -6;1 -5 3];b=[-1;-2;1];>> rank(A), rank([A,b])ans =3ans =3可见方程组唯一解>> x=A\bx =-0.4706-0.2941(3)>> A=[4 1;3 2;1 -5];b=[1;1;1];>> rank(A), rank([A,b])ans =2ans =3可见方程组无解>> x=A\bx =0.3311-0.1219最小二乘近似解(4)>> a=[2,1,-1,1;1,2,1,-1;1,1,2,1];b=[1 2 3]';%注意b的写法>> rank(a),rank([a,b])ans =3ans =3 rank(a)==rank([a,b])<4说明有无穷多解>> a\bans =110一个特解Page65 Ex3>> a=[2,1,-1,1;1,2,1,-1;1,1,2,1];b=[1,2,3]'; >> x=null(a),x0=a\bx =-0.62550.6255-0.20850.4170x0 =11通解kx+x0Page65 Ex 4>> x0=[0.2 0.8]';a=[0.99 0.05;0.01 0.95]; >> x1=a*x, x2=a^2*x, x10=a^10*x>> x=x0;for i=1:1000,x=a*x;end,xx =0.83330.1667>> x0=[0.8 0.2]';>> x=x0;for i=1:1000,x=a*x;end,xx =0.83330.1667>> [v,e]=eig(a)v =0.9806 -0.70710.1961 0.7071e =1.0000 00 0.9400>> v(:,1)./xans =1.17671.1767成比例,说明x是最大特征值对应的特征向量Page65Ex5用到公式(3.11)(3.12)>> B=[6,2,1;2.25,1,0.2;3,0.2,1.8];x=[25 5 20]';>> C=B/diag(x)C =0.2400 0.4000 0.05000.0900 0.2000 0.01000.1200 0.0400 0.0900>> A=eye(3,3)-CA =0.7600 -0.4000 -0.0500-0.0900 0.8000 -0.0100-0.1200 -0.0400 0.9100>> D=[17 17 17]';x=A\Dx =37.569625.786224.7690Page65Ex 6(1)>> a=[4 1 -1;3 2 -6;1 -5 3];det(a),inv(a),[v,d]=eig(a) ans =-94ans =0.2553 -0.0213 0.04260.1596 -0.1383 -0.22340.1809 -0.2234 -0.0532v =0.0185 -0.9009 -0.3066-0.7693 -0.1240 -0.7248-0.6386 -0.4158 0.6170d =-3.0527 0 00 3.6760 00 0 8.3766(2)>> a=[1 1 -1;0 2 -1;-1 2 0];det(a),inv(a),[v,d]=eig(a) ans =1ans =2.0000 -2.0000 1.00001.0000 -1.0000 1.00002.0000 -3.0000 2.0000v =-0.5773 0.5774 + 0.0000i 0.5774 - 0.0000i -0.5773 0.5774 0.5774-0.5774 0.5773 - 0.0000i 0.5773 + 0.0000id =1.0000 0 00 1.0000 + 0.0000i 00 0 1.0000 - 0.0000i(3)>> A=[5 7 6 5;7 10 8 7;6 8 10 9;5 7 9 10]A =5 76 57 10 8 76 8 10 95 7 9 10>> det(A),inv(A), [v,d]=eig(A)ans =1ans =68.0000 -41.0000 -17.0000 10.0000 -41.0000 25.0000 10.0000 -6.0000 -17.0000 10.0000 5.0000 -3.0000 10.0000 -6.0000 -3.0000 2.0000v =0.8304 0.0933 0.3963 0.3803-0.5016 -0.3017 0.6149 0.5286-0.2086 0.7603 -0.2716 0.55200.1237 -0.5676 -0.6254 0.5209d =0.0102 0 0 00 0.8431 0 00 0 3.8581 00 0 0 30.2887(4)(以n=5为例)方法一(三个for)n=5;for i=1:n, a(i,i)=5;endfor i=1:(n-1),a(i,i+1)=6;endfor i=1:(n-1),a(i+1,i)=1;enda方法二(一个for)n=5;a=zeros(n,n);a(1,1:2)=[5 6];for i=2:(n-1),a(i,[i-1,i,i+1])=[1 5 6];enda(n,[n-1 n])=[1 5];a方法三(不用for)n=5;a=diag(5*ones(n,1));b=diag(6*ones(n-1,1));c=diag(ones(n-1,1));a=a+[zeros(n-1,1),b;zeros(1,n)]+[zeros(1,n);c,zeros(n-1,1)]下列计算>> det(a)ans =665>> inv(a)ans =0.3173 -0.5865 1.0286 -1.6241 1.9489-0.0977 0.4887 -0.8571 1.3534 -1.62410.0286 -0.1429 0.5429 -0.8571 1.0286-0.0075 0.0376 -0.1429 0.4887 -0.58650.0015 -0.0075 0.0286 -0.0977 0.3173 >> [v,d]=eig(a)v =-0.7843 -0.7843 -0.9237 0.9860 -0.9237 0.5546 -0.5546 -0.3771 -0.0000 0.3771 -0.2614 -0.2614 0.0000 -0.1643 0.0000 0.0924 -0.0924 0.0628 -0.0000 -0.0628 -0.0218 -0.0218 0.0257 0.0274 0.0257d =0.7574 0 0 0 00 9.2426 0 0 00 0 7.4495 0 00 0 0 5.0000 00 0 0 0 2.5505Page65 Ex 7(1)>> a=[4 1 -1;3 2 -6;1 -5 3];[v,d]=eig(a) v =0.0185 -0.9009 -0.3066-0.7693 -0.1240 -0.7248-0.6386 -0.4158 0.6170d =-3.0527 0 00 3.6760 00 0 8.3766>> det(v)ans =-0.9255 %v行列式正常,特征向量线性相关,可对角化>> inv(v)*a*v验算ans =-3.0527 0.0000 -0.00000.0000 3.6760 -0.0000-0.0000 -0.0000 8.3766>> [v2,d2]=jordan(a)v2 =0.0798 0.0076 0.91270.1886 -0.3141 0.1256-0.1605 -0.2607 0.4213d2 =8.3766 0 00 -3.0527 - 0.0000i 00 0 3.6760 + 0.0000i>> v2\a*v2ans =8.3766 0 0.00000.0000 -3.0527 0.00000.0000 0.0000 3.6760也可用jordan特征向量不同>> v(:,1)./v2(:,2)对应相同特征值的特征向量成比例ans =2.44912.44912.4491(2)>> a=[1 1 -1;0 2 -1;-1 2 0];[v,d]=eig(a)v =-0.5773 0.5774 + 0.0000i 0.5774 - 0.0000i-0.5773 0.5774 0.5774-0.5774 0.5773 - 0.0000i 0.5773 + 0.0000id =1.0000 0 00 1.0000 + 0.0000i 00 0 1.0000 - 0.0000i>> det(v)ans =-5.0566e-028-5.1918e-017i v的行列式接近0,角化>> [v,d]=jordan(a)v =1 0 1特征向量线性相关,不可对1 0 01 -1 0d =1 1 00 1 10 0 1 jordan标准形不是对角的,所以不可对角化(3)>> A=[5 7 6 5;7 10 8 7;6 8 10 9;5 7 9 10]A =5 76 57 10 8 76 8 10 95 7 9 10>> [v,d]=eig(A)v =0.8304 0.0933 0.3963 0.3803-0.5016 -0.3017 0.6149 0.5286-0.2086 0.7603 -0.2716 0.55200.1237 -0.5676 -0.6254 0.5209d =0.0102 0 0 00 0.8431 0 00 0 3.8581 00 0 0 30.2887>> inv(v)*A*vans =0.0102 0.0000 -0.0000 0.00000.0000 0.8431 -0.0000 -0.0000-0.0000 0.0000 3.8581 -0.0000-0.0000 -0.0000 0 30.2887本题用jordan不行,原因未知(4)参考6(4)和7(1)Page65Exercise 8只有(3)对称,且特征值全部大于零,所以是正定矩阵. Page65Exercise 9(1)>> a=[4 -3 1 3;2 -1 3 5;1 -1 -1 -1;3 -2 3 4;7 -6 -7 0] >> rank(a)ans =3>> rank(a(1:3,:))ans =2>> rank(a([1 2 4],:)) 1,2,4行为最大无关组ans =3>> b=a([1 2 4],:)';c=a([3 5],:)';>> b\c线性表示的系数ans =0.5000 5.0000-0.5000 1.00000 -5.0000Page65Exercise 10>> a=[1 -2 2;-2 -2 4;2 4 -2]>> [v,d]=eig(a)v =0.3333 0.9339 -0.12930.6667 -0.3304 -0.6681-0.6667 0.1365 -0.7327d =-7.0000 0 00 2.0000 00 0 2.0000>> v'*vans =1.0000 0.0000 0.00000.0000 1.0000 00.0000 0 1.0000 v确实是正交矩阵Page65 Exercise 11设经过6个电阻的电流分别为i1, ..., i6.列方程组如下20-2i1=a; 5-3i2=c; a-3i3=c; a-4i4=b; c-5i5=b; b-3i6=0;i1=i3+i4;i5=i2+i3;i6=i4+i5;计算如下>> A=[1 0 0 2 0 0 0 0 0;0 0 1 0 3 0 0 0 0;1 0 -1 0 0 -3 0 0 0; 1 -1 0 00 0 -4 0 0;0 -1 1 0 0 0 0 -5 0;0 1 0 0 0 0 0 0 -3; 0 0 0 1 0 -1 -1 0 0;0 0 0 0 -1 -10 1 0;0 0 0 0 0 0 -1 -1 1];>>b=[20 5 0 0 0 0 0 0 0]'; A\bans =13.34536.44018.54203.3274-1.18071.60111.72630.42042.1467Page65Exercise 12>> A=[1 2 3;4 5 6;7 8 0];>> left=sum(eig(A)), right=sum(trace(A))left =6.0000right =6>> left=prod(eig(A)), right=det(A)原题有错, (-1)^n应删去left =27.0000right =27>> fA=(A-p(1)*eye(3,3))*(A-p(2)*eye(3,3))*(A-p(3)*eye(3,3))fA = 1.0e-012 *0.0853 0.1421 0.02840.1421 0.1421 0-0.0568 -0.1137 0.1705>> norm(fA) f(A)范数接近0ans =2.9536e-013Chapter 4Page84 Exercise 1(1)roots([1 1 1])(2)roots([3 0 -4 0 2 -1])(3)p=zeros(1,24);p([1 17 18 22])=[5 -6 8 -5];roots(p)(4)p1=[2 3];p2=conv(p1, p1);p3=conv(p1, p2);p3(end)=p3(end)-4; %原p3最后一个分量-4roots(p3)Page84 Exercise 2fun=inline('x*log(sqrt(x^2-1)+x)-sqrt(x^2-1)-0.5*x');fzero(fun,2)Page84Exercise 3fun=inline('x^4-2^x');fplot(fun,[-2 2]);grid on;fzero(fun,-1),fzero(fun,1),fminbnd(fun,0.5,1.5)Page84 Exercise 4fun=inline('x*sin(1/x)','x');fplot(fun, [-0.1 0.1]);x=zeros(1,10);for i=1:10, x(i)=fzero(fun,(i-0.5)*0.01);end;x=[x,-x]Page84 Exercise 5fun=inline('[9*x(1)^2+36*x(2)^2+4*x(3)^2-36;x(1)^2-2*x(2)^2-20*x(3);16*x(1)-x(1)^3-2*x(2)^2-16*x(3)^2]','x');[a,b,c]=fsolve(fun,[0 0 0])Page84Exercise 6fun=@(x)[x(1)-0.7*sin(x(1))-0.2*cos(x(2)),x(2)-0.7*cos(x(1))+0.2*sin(x(2))];[a,b,c]=fsolve(fun,[0.5 0.5])Page84Exercise 7clear; close; t=0:pi/100:2*pi;x1=2+sqrt(5)*cos(t); y1=3-2*x1+sqrt(5)*sin(t);x2=3+sqrt(2)*cos(t); y2=6*sin(t);plot(x1,y1,x2,y2); grid on;作图发现4个解的大致位置,然后分别求解y1=fsolve('[(x(1)-2)^2+(x(2)-3+2*x(1))^2-5,2*(x(1)-3)^2+(x(2)/3)^2-4]',[1.5,2]) y2=fsolve('[(x(1)-2)^2+(x(2)-3+2*x(1))^2-5,2*(x(1)-3)^2+(x(2)/3)^2-4]',[1.8,-2]) y3=fsolve('[(x(1)-2)^2+(x(2)-3+2*x(1))^2-5,2*(x(1)-3)^2+(x(2)/3)^2-4]',[3.5,-5]) y4=fsolve('[(x(1)-2)^2+(x(2)-3+2*x(1))^2-5,2*(x(1)-3)^2+(x(2)/3)^2-4]',[4,-4]) Page84Exercise 8(1)clear;fun=inline('x.^2.*sin(x.^2-x-2)');fplot(fun,[-2 2]);grid on;作图观察x(1)=-2;x(3)=fminbnd(fun,-1,-0.5);x(5)=fminbnd(fun,1,2);fun2=inline('-x.^2.*sin(x.^2-x-2)');x(2)=fminbnd(fun2,-2,-1);x(4)=fminbnd(fun2,-0.5,0.5);x(6)=2feval(fun,x)答案:以上x(1)(3)(5)是局部极小,x(2)(4)(6)是局部极大,从最后一句知道x(1)全局最小,x(2)最大。