数学建模实验小结由刀豆文库小编整理,希望给你工作、学习、生活带来方便,猜你可能喜欢“数学建模实验总结”。
例1-1 >> r=2;V=4/3*pi*r^3 V =
33.5103 例2-1 计算s=...>> s=0;>> for n=1:100 s=s+1/n/n;end >> s s =
1.6350 例2-5 两个一元函数y=x3-x-1,y=x .2sin(5x)在区间-1
y=abs(x).^0.2.*sin(5*x);plot(x,y,':ro');hold off;
曲面图 >> xa=6:8;ya=1:4;>> [x,y]=meshgrid(xa,ya);>> z=x.^2+y.^2;>> mesh(x,y,z)>> [x,y,z] ans =
例2-6 二元函数图z=xexp(-x2-y2).xa=-2:0.2:2;ya=xa;
[x,y]=meshgrid(xa,ya);z=x.*exp(-x.^2-y.^2);mesh(x,y,z);pause;surf(x,y,z);pause;
contour(x,y,z,[0.1,0.1]);pause mesh(x,y,z);
Page40 1.先在编辑器窗口写下列M函数,保存为ex2_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.4000 s =
12.1124 2.s=log(1);n=0;while s
s=s+log(1+n);end
m=n 3.F(1)=1;F(2)=1;k=2;x=0;e=1e-8;a=(1+sqrt(5))/2;while abs(x-a)>e
k=k+1;F(k)=F(k-1)+F(k-2);x=F(k)/F(k-1);end
a,x,k m =
a =
1.6180 x =
1.6180 k = 4.clear;tic;s=0;for i=1:1000000 s=s+sqrt(3)/2^i;end
s,toc
tic;s=0;i=1;
while i
s=s+sqrt(3)/2^i;i=i+1;
end
s,toc tic;s=0;
i=1:1000000;
s=sqrt(3)*sum(1./2.^i);s,toc
s =
1.7321 Elapsed time is 2.038973 seconds.s =
1.7321 Elapsed time is 2.948968 seconds.s =
1.7321 Elapsed time is 0.453414 seconds 5.t=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)
6.(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)
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)%page41, ex7 x=-1.5:0.05:1.5;
y=1.1*(x>1.1)+x.*(x=-1.1)-1.1*(x
Page59 1.>> a=[1 2 3];b=[2 4 3];>> a./b ans =
0.5000
0.5000
1.0000 >> a.b ans = >> a/b ans =
0.6552 >> ab ans =
0
0
0
0
0
0
0.6667
1.3333
1.0000 2.(1)>> a=[4 1-1;3 2-6;1-5 3];b=[9;-2;1];>> ab ans =
2.3830
1.4894 2.0213(2)>> a=[4-3 3;3 2-6;1-5 3],b=[-1;-2;1] a =
-5b =
>> ab ans =
-0.4706
-0.2941
0(3)>> a=[4 1;3 2;1-5],b=[1;1;1] a =
-5 b =
>> ab ans =
0.3311
-0.1219(4)>> a=[2 1-1 1;1 2 1-1;1 1 2 1],b=[1;2;3] a =
-1b =
>> ab ans =
0
0 6.(1)>> a=[4 1-1;3 2-6;1-5 3];>> b=det(a),inv(a),[V,D]=eig(a)b =
-94 ans =
0.2553
-0.0213
0.0426
0.1596
-0.1383
-0.2234
0.1809
-0.2234
-0.0532 V =
0.0185
-0.9009
-0.3066
-0.7693
-0.1240
-0.7248
-0.6386
-0.4158
0.6170 D =
-3.0527
0
0
0
3.6760
0
0
0
8.3766(2)>> a=[1 1-1;0 2-1;-1 2 0];b=det(a),inv(a),[V,D]=eig(a)b =
ans =
2.0000
-2.0000
1.0000
1.0000
-1.0000
1.0000
2.0000
-3.0000
2.0000 V =
-0.5773
0.5774 + 0.0000i
0.57740.0000i
0.5773 + 0.0000i D =
1.0000
0
0
0
1.0000 + 0.0000i
0
0
0
1.00000.0000i
-0.5773
0.5774
0.5774
-0.5774
0.57730.0000i >> det(V)ans =
-5.0566e-028-5.1918e-017i
%V的行列式接近0, 特征向量线性相关,不可对角化(3)>> a=[5 7 6 5;7 10 8 7;6 8 10 9;5 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.5520
0.1237
-0.5676
-0.6254
0.5209 D =
0.0102
0
0
0
0
0.8431
0
0
0
0
3.8581
0
0
0
0
30.2887 >> inv(V)*a*V ans =
0.0102
0.0000
-0.0000
0.0000
0.0000
0.8431
-0.0000
-0.0000
-0.0000
0.0000
3.8581
-0.0000
-0.0000
-0.0000
0
30.2887 8
对称阵A为正定的充分必要条件是:A的特征值全为正。只有(3)对称, 且特征值全部大于零, 所以(3)是正定矩阵.例4.2用2次多项式拟合下列数据。>> clear;x=[0.1,0.2,0.15,0,-0.2,0.3];>> y=[0.95,0.84,0.86,1.06,1.50,0.72];>> p=polyfit(x,y,2)p =
1.7432
-1.6959
1.0850 得到二次拟合式:1.7432x^2-1.6959x+1.0850 >> xi=-0.2:0.01:0.3;>> yi=polyval(p,xi);plot(x,y,'o',xi,yi);
例4.3 求函数y=x*sin(x^2-x-1)在(-2,-0.1)内的零点。>> fun=inline('x*sin(x^2-x-1)','x')fun =
Inline function:
fun(x)= x*sin(x^2-x-1)>> fzero(fun,[-2,-0.1])??? Error using ==> fzero at 292 The function values at the interval endpoints must differ in sign.>> fplot(fun,[-2,-0.1]);grid on;
>> [x,f,h]=fsolve(fun,-1.6),[x,f,h]=fsolve(fun,-0.6)Optimization terminated: first-order optimality is le than options.TolFun.x =
-1.5956 f =
1.4909e-009 h =
Optimization terminated: first-order optimality is le than options.TolFun.x =
-0.6180 f =
-3.3152e-012 h =
例4.4求下列方程组在原点附近的解
>> fun=inline('[4*x(1)-x(2)+exp(x(1))/10-1,-x(1)+4*x(2)+x(1)^2/8]','x');[x,f,h]=fsolve(fun,[0,0])Optimization terminated: first-order optimality is le than options.TolFun.x =
0.2326
0.0565 f =
1.0e-006 *
0.0908
0.1798 h =
例4.5 求二元函数f(x,y)=5-x^4-y^4+4*x*y在原点附近的极大值。(等价于求-f(x,y)的极小值)
>> fun=inline('x(1)^4+x(2)^4-4*x(1)*x(2)-5');>> [x,g]=fminsearch(fun,[0,0])x =
1.0000
1.0000 g =
-7.0000 例4.6 用Newton迭代法求下列方程的正根,要求精度为10的-6次 X^2-3x+e^x=2 >> fun=inline('x^2-3*x+exp(x)-2');>> fplot(fun,[0,2]);>> grid on;
%M函数 newton.m function x =newton(fname,dfname,x0,e)if nargine
x0=x;x=x0-feval(fname,x0)/feval(dfname,x0);end
>> dfun=inline('2*x-3+exp(x)');format long;newton(fun,dfun,1.5,1e-6),format short ans =
1.*** 例4.7 用函数y=a*e^(b*x)拟合例4.2的数据。>> fun=inline('c(1)*exp(c(2)*x)','c','x');>> x=[0.1,0.2,0.15,0,-0.2,0.3];y=[0.95,0.84,0.86,1.06,1.50,0.72];>> c=lsqcurvefit(fun,[0,0],x,y)Optimization terminated: relative function value changing by le than OPTIONS.TolFun.c =
1.0997
-1.4923
PAGE 77 1.%Exercise 1(1)roots([1 1 1])%Exercise 1(2)
roots([3 0-4 0 2-1])%Exercise 1(3)p=zeros(1,24);
p([1 17 18 22])=[5-6 8-5];roots(p)
%Exercise 1(4)p1=[2 3];
p2=conv(p1, p1);p3=conv(p1, p2);
p3(end)=p3(end)-4;%原p3最后一个分量-4 roots(p3)2.%Exercise 2
fun=inline('x*log(sqrt(x^2-1)+x)-sqrt(x^2-1)-0.5*x');fzero(fun,2)3.%Exercise 3
fun=inline('x^4-2^x');fplot(fun,[-2 2]);grid on;
fzero(fun,-1),fzero(fun,1),fminbnd(fun,0.5,1.5)
4.%Exercise 4
fun=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] x =
Columns 1 through 11
0.0050
0.0152
0.0245
0.0354
0.0455
0.0531
0.0637
0.0796
0.0796
0.1061
-0.0050
Columns 12 through 20
-0.0152
-0.0245
-0.0354
-0.0455
-0.0531
-0.0637
-0.0796
-0.0796
-0.1061
5.%Exercise 5
fun=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])6.%Exercise 6
fun=@(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])7.%Exercise 7
clear;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])
8.%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)=2 feval(fun,x)x =
-2.0000
-1.5326
-0.7315
-0.0000
1.5951
2.0000 ans =
-3.0272
2.2364
-0.3582
-0.0000
-2.2080
0
%答案: 以上x(1)(3)(5)是局部极小,x(2)(4)(6)是局部极大,从最后一句知道x(1)全局最小,x(2)最大。
%Exercise 8(2)clear;
fun=inline('3*x.^5-20*x.^3+10');fplot(fun,[-3 3]);grid on;%作图观察
x(1)=-3;
x(3)=fminsearch(fun,2.5);
fun2=inline('-(3*x.^5-20*x.^3+10)');x(2)=fminsearch(fun2,-2.5);x(4)=3;feval(fun,x)ans =
-179
-54
199
%Exercise 8(3)
fun=inline('abs(x^3-x^2-x-2)');fplot(fun,[0 3]);grid on;%作图观察
fminbnd(fun,1.5,2.5)
fun2=inline('-abs(x^3-x^2-x-2)');fminbnd(fun2,0.5,1.5)ans =
2.0000 ans =
1.0000
9.%Exercise 9 close;
x=-2:0.1:1;y=-7:0.1:1;[x,y]=meshgrid(x,y);
z=y.^3/9+3*x.^2.*y+9*x.^2+y.^2+x.*y+9;mesh(x,y,z);grid on;%作图观察
fun=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(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])%求极大值
x =
0
0 x =
-0.3333
-6.0000
10.clear;t=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];p2=polyfit(t,c,2)p3=polyfit(t,c,3)
fun=inline('a(1)*exp(a(2)*(t-14).^2)','a','t');
a=lsqcurvefit(fun,[0 0],t,c)%初值可以试探 f=feval(fun, a,t)
norm(f-c)%拟合效果
plot(t,c,t,f)%作图检验
fun2=inline('b(1)*sin(pi/12*t+b(2))+20','b','t');%原题修改f(x)+20 b=lsqcurvefit(fun2,[0 0],t,c)figure
f2=feval(fun2, b,t)
norm(f2-c)%拟合效果
plot(t,c,t,f2)%作图检验
Page 94 chapter5 1.x=[0 4 10 12 15 22 28 34 40];y=[0 1 3 6 8 9 5 3 0];trapz(x,y)ans =
178.5000 2.>> x=[0 4 10 12 15 22 28 34 40];y=[0 1 3 6 8 9 5 3 0];diff(y)./diff(x)
ans =
0.2500
0.3333
1.5000
0.6667
0.1429
-0.6667
-0.3333
-0.5000 3.xa=-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 4.t=0:0.01:1.5;x=log(cos(t));y=cos(t)-t.*sin(t);dydx=gradient(y,x)
[x_1,id]=min(abs(x-(-1)));%找最接近x=-1的点 dydx(id)5.(1)(2)fun=inline('exp(2*x).*cos(x).^3');quadl(fun,0,2*pi)(3)fun=@(x)x.*log(x.^4).*asin(1./x.^2);quadl(fun,1,3)(4)fun=@(x)sin(x)./x;
quadl(fun,1e-10,1)%注意由于下限为0,被积函数没有意义,用很小的1e-10代替(5)(6)fun=inline('sqrt(1+r.^2.*sin(th))','r','th');dblquad(fun,0,1,0,2*pi)(7)首先建立84页函数dblquad2 clear;
fun=@(x,y)1+x+y.^2;clo=@(x)-sqrt(2*x-x.^2);dup=@(x)sqrt(2*x-x.^2);dblquad2(fun,0,2,clo,dhi,100)%Exercise 6
t=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)10(1)(2)
%先在程序编辑器,写下列函数,保存为ex5_10_2f function d=ex5_10_2f(fname,a,h0,e)
h=h0;d=(feval(fname,a+h)-2*feval(fname,a)+feval(fname,a-h))/(h*h);d0=d+2*e;
while abs(d-d0)>e d0=d;h0=h;h=h0/2;
d=(feval(fname,a+h)-2*feval(fname,a)+feval(fname,a-h))/(h*h);end %再在指令窗口执行
fun=inline('x.^2*sin(x.^2-x-2)','x');d=ex5_10_2f(fun,1.4,0.1,1e-3)13.fun=inline('5400*v./(8.276*v.^2+2000)','v');quadl(fun,15,30)