Matlab软件与基础数学实验程序部分.docx
- 文档编号:28750863
- 上传时间:2023-07-19
- 格式:DOCX
- 页数:23
- 大小:952.53KB
Matlab软件与基础数学实验程序部分.docx
《Matlab软件与基础数学实验程序部分.docx》由会员分享,可在线阅读,更多相关《Matlab软件与基础数学实验程序部分.docx(23页珍藏版)》请在冰豆网上搜索。
Matlab软件与基础数学实验程序部分
追击问题:
一敌舰在某海域内以椭圆轨迹航行,其在时间t时刻的坐标为:
x(t)=10+20cost
y(t)=20+5sint
我方战舰恰位于原点处,我战舰向敌舰发射制导鱼雷,鱼雷的速率为20,其运行方向始终指向敌舰,试问敌舰航行在何处将被击中?
2.若敌舰的运行轨迹变为
x(t)=10+20cost
y(t)=20+20sint
试问敌舰航行在何处将被击中?
(无法击中)
3.若敌舰的运行轨迹变为
x(t)=10+20cost
y(t)=20+20sint
鱼雷速率提高至21,结果如何?
%Matlab程序:
clear;clc
h=;%时间步长
k=1;
t
(1)=0;x
(1)=0;y
(1)=0;%初始值
r=10;
whiler>=%k<=250%
m=(20+5*sin(t(k))-y(k))/(10+20*cos(t(k))-x(k)++;
if10+20*cos(t(k))-x(k)>=0
x(k+1)=x(k)+20*h/sqrt(1+m^2);
else
x(k+1)=x(k)-20*h/sqrt(1+m^2);
end
if20+5*sin(t(k))-y(k)>=0
y(k+1)=y(k)+20*h/sqrt(1+1/m/m);
else
y(k+1)=y(k)-20*h/sqrt(1+1/m/m);
end
r=(x(k)-10-20*cos(t(k)))^2+(y(k)-20-5*sin(t(k)))^2;
r=sqrt(r);
t(k+1)=h*k;
k=k+1;
plot(10+20*cos(t(k)),20+5*sin(t(k)),'r*')
holdon
axis([-1032-330]);
plot(x,y,'o')
pause
end
t=t(end),x=x(end),y=y(end)
t=
x=
y=
第二问:
速度相同无法击中
第三问:
t=
x=
y=
%Matlab程序:
clear;clc
h=;%时间步长
k=1;
t
(1)=0;x
(1)=0;y
(1)=0;%初始值
r=10;
whiler>=%k<=250%
m=(20+20*sin(t(k))-y(k))/(10+20*cos(t(k))-x(k)++;
if10+20*cos(t(k))-x(k)>=0
x(k+1)=x(k)+22*h/sqrt(1+m^2);
else
x(k+1)=x(k)-22*h/sqrt(1+m^2);
end
if20+20*sin(t(k))-y(k)>=0
y(k+1)=y(k)+22*h/sqrt(1+1/m/m);
else
y(k+1)=y(k)-22*h/sqrt(1+1/m/m);
end
r=(x(k)-10-20*cos(t(k)))^2+(y(k)-20-20*sin(t(k)))^2;
r=sqrt(r);
t(k+1)=h*k;
k=k+1;
plot(10+20*cos(t(k)),20+20*sin(t(k)),'r*')
holdon
axis([-1232-242]);
plot(x,y,'o')
pause
end
t=t(end),x=x(end),y=y(end)
课本P81
1.某农夫有一个半径10米的圆形牛栏,长满了草.他要将一头牛栓在牛栏边界的栏桩上,但只让牛吃到一半草,问栓牛鼻的绳子应为多长?
设拴牛的绳子长为r,以圆形牛栏C1的圆心为原点建立直角坐标系,见图1,不妨设拴牛的栏桩为图1中圆形牛栏C1上的B点,其坐标为(10,0),则所求问题转化为:
求出r,使得以B点为圆心,半径为r的圆C2与圆C1相交部分的面积
是圆C1面积的一半。
解法一:
由于圆形牛栏C1和圆C2的方程分别为:
C1:
x2+y2=100
C2:
(x-10)2+y2=r2
(1)
联立方程C1,C2,可得两交点分别为:
,
设牛吃草的面积为S,即圆C1与C2的相交部分,则根据题意,S应为圆C1面积的一半,即
由图可知,S的面积可由下面的定积分计算得到:
对于上式中的积分
,令u=x-10,则上式可化简为
S
上式通过简单积分运算可化为
f=inline('-r^3/800*sqrt(400-r^2)*r^2*asin(r/20)+*pi*r^*sqrt(400*r^2-r^4)+r^2/800*sqrt(400*r^2-r^4)-50*asin(1-r^2/200)','r');
fzero(f,10)
ans=
另解:
S0=1/2*pi*10^2;
r=10;%给定r的迭代初值;
S=314;
whileabs(S-S0)>
r=r+;%r的值增加;
u=-r:
:
-r^2/20;
s1=trapz(u,sqrt(r^2-u.^2));
x=10-r^2/20:
:
10;
s2=trapz(x,sqrt(10^2-x.^2));
S=2*(s1+s2);
end%当误差小于eps时,循环结束;
error=abs(S0-S)%显示误差;
[r,S]%显示r的值和面积;
error=
ans=
11.5900
解法二(几何方法):
x=10:
:
14;
fx=*x.*x.*acos*x)+50*acos*x.*x)*x.*sqrt*x.*x)-25*pi;
plot(x,fx)
grid
y=inline('*x.*x.*acos(x/20)+50*acos(1-x.*x/200)*x.*sqrt*x.*x)-25*pi')
fzero(y,[10,14])
ans=
2.如图所示,为了在海岛I与某城市C之间铺设一条地下光缆,每千米光缆铺设成本在水下部分是C1,在地下的部分是C2,为使的铺设该光缆的总成本最低,光缆C1的转折点P(海岸线上)应取在何处?
如果实际测得海岛I与城市C之间水平距离L=30km,海岛距海岸垂直距离h1=15km,城市距海岸线垂直距离h=10km,C1=3000万元/km,C2=1500万元/km,求p点坐标(误差<10-3km)
ezplot('3000*sqrt((30-x)^2+225)+1500*sqrt(x^2+100)',[0,30])
symsx
z=3000*sqrt((30-x)^2+225)+1500*sqrt(x^2+100);
f=inline('1500/(1125-60*x+x^2)^(1/2)*(-60+2*x)+1500/(x^2+100)^(1/2)*x');
a=20;
b=25;
dlt=;
k=1;
whileabs(b-a)>dlt
c=(a+b)/2;
iff(c)==0
break;
elseiff(c)*f(b)<0
a=c;
elseb=c;
end
k=k+1
end%二分法求根
x=c%距城市C的水平位置;
L=eval(z)%最低成本
x=
L=+004
3.有一艘宽为5m的长方形驳船,欲过某河道的直角弯,经测量知河道的宽度10m和12m,如图所示,设问,驳船要驶过直角弯,驳船的长度不能超过多少米?
(误差<
m)
ezplot('10/sin(x)+12/cos(x)-5*tan(x)-5/tan(x)',[0,pi/2])
symsx
z=10/sin(x)+12/cos(x)-5*tan(x)-5/tan(x);
f=inline('-10/sin(x)^2*cos(x)+12/cos(x)^2*sin(x)-5-5*tan(x)^2+5/tan(x)^2*(1+tan(x)^2)');
a=;
b=;
dlt=;
k=1;
whileabs(b-a)>dlt
c=(a+b)/2;
iff(c)==0
break;
elseiff(c)*f(b)<0
a=c;
elseb=c;
end
k=k+1
end%二分法求根
x=c%角度值
L=eval(z)%驳船的长度
输出结果:
x=
L=
4.一个对称的地下油库,横截面为圆,中心位置上的半径为3m,上下底上的半径为2m,高为12m,纵截面的两侧是顶点在中心位置的抛物线。
(1)试求:
油库内油面的深度为h(从底部算起)时,库内油量的容积V(h);
(2)根据刻度能读出油库内油量的多少,试给出V=20,30的高度。
建立坐标系求出抛物线方程:
抛物线方程可设为:
x=ay2+by+c,通过A(3,0)、B(2,6)、C(2,-6)
得:
程序
clc;clear;
h=input('Pleaseinputh=');
ifh>=0&h<=12%转化为标准高度!
h=h-6;
v=pi*(h^5/6480-h^3/18+9*h+216/5);
fprintf('v=%.4f\n',v);
else
fprintf('输入数据超出高度范围');
end
(2)由体积反求高度
clear;clc;
f20=inline('pi*(h^5/6480-h^3/18+9*h+216/5)-20','h');
f30=inline('pi*(h^5/6480-h^3/18+9*h+216/5)-30','h');
h20=fzero(f20,[-6,6])+6
h30=fzero(f30,[-6,6])+6
h20=
h30=
课本P91
1.Feigenbaum在做研究时,对超越函数
(
为非负实数)进行了分岔与混沌的研究,试利用迭代格式
,做出相应的Feigenbaum图。
clear;clc;
holdon
axis([0,4,-4,4]);
holdon
forr=0:
:
x=[];
fork=2:
150
x(k)=r*sin(pi*x(k-1));
end
%pause
fork=101:
105
plot(r,x(k),'k.');
end
end
第二题
2.(Henon吸引子)混沌和分形的著名例子,迭代模型为
取初值x0=0,y0=0,进行3000次迭代,对于k>1000,在(xk,yk)处画一点(注意不要连线)可得所谓Henon引力线图.试写出程序,画出图形.
clear;clc;
holdon
x=[0];
y=[0];
k=1;
whilek<=3000
a=x(k);
x(k+1)=1+y(k)*x(k)*x(k);
y(k+1)=*a;
k=k+1;
ifk>1000
plot(x(k),y(k),'b.','Markersize',5)
end
end
课本P101
4.椭圆周长计算
梯形公式计算:
clear;
f=inline('sqrt(1+3.*sin(t).^2)');
a=0;b=2.*pi;n=1;
h=(b-a)/n;
t1=h/2*(f(a)+f(b));
er=1;k=1;
whileer>
s=0;
fori=1:
n
s=s+f(a+(i-1/2)*h);
end
t2=(t1+h*s)/2;
er=abs(t2-t1);
fprintf('n=%.0f,p=%.6f,r=%.6f\n',k,t2,er);
n=2*n;h=h/2;t1=t2;
k=k+1;
end
n=1,p=,r=
这是为什么?
?
?
clear;
f=inline('sqrt(1+3.*sin(t).^2)');
a=0;b=pi;n=1;
h=(b-a)/n;
t1=h/2*(f(a)+f(b));
er=1;k=1;
whileer>
s=0;
fori=1:
n
s=s+f(a+(i-1/2)*h);
end
t2=(t1+h*s)/2;
er=abs(t2-t1);
fprintf('n=%.0f,p=%.6f,r=%.6f\n',k,t2,er);
n=2*n;h=h/2;t1=t2;
k=k+1;
end
t2=t2*2
n=1,p=,r=
n=2,p=,r=
n=3,p=,r=
n=4,p=,r=
n=5,p=,r=
t2=
clear;
f=inline('sqrt(1+.*x.*x./(4-x.*x))');
a=0;b=2.;n=1;
h=(b-a)/n;
t1=h/2*(f(a)+f(b));
er=1;k=1;
whileer>
s=0;
fori=1:
n
s=s+f(a+(i-1/2)*h);
end
t2=(t1+h*s)/2;
er=abs(t2-t1);
fprintf('n=%.0f,p=%.6f,r=%.6f\n',k,t2,er);
n=2*n;h=h/2;t1=t2;
k=k+1;
end
t2=t2*4
Warning:
Dividebyzero.
(Type"warningoffMATLAB:
divideByZero"tosuppressthiswarning.)
>InC:
\MATLAB6p5\toolbox\matlab\funfun\atline13
InC:
\MATLAB6p5\toolbox\matlab\funfun\@inline\atline25
n=1,p=Inf,r=NaN
t2=
Inf
clear;
f=inline('sqrt(1+.*x.*x./(4-x.*x))');
a=0;b=1.;n=1;
h=(b-a)/n;
t1=h/2*(f(a)+f(b));
er=1;k=1;
whileer>
s=0;
fori=1:
n
s=s+f(a+(i-1/2)*h);
end
t2=(t1+h*s)/2;
er=abs(t2-t1);
fprintf('n=%.0f,p=%.6f,r=%.6f\n',k,t2,er);
n=2*n;h=h/2;t1=t2;
k=k+1;
end
t2=t2*4
n=1,p=,r=
n=2,p=,r=
n=3,p=,r=
n=4,p=,r=
n=5,p=,r=
n=6,p=,r=
n=7,p=,r=
n=8,p=,r=
n=9,p=,r=
n=10,p=,r=
n=11,p=,r=
n=12,p=,r=
n=13,p=,r=
n=14,p=,r=
n=15,p=,r=
n=16,p=,r=
n=17,p=,r=
n=18,p=,r=
n=19,p=,r=
n=20,p=,r=
n=21,p=,r=
n=22,p=,r=
n=23,p=,r=
.....
ans=
注意:
计算量太大!
方法二
clear;
f=inline('sqrt(1+.*x.*x./(4-x.*x))');
a=0;b=;
z=4*quad(f,a,b,;
fprintf('z=%.4f\n',z)
或:
clear;
f=inline('sqrt(1+3.*sin(t).^2)');
a=0;b=2.*pi;n=1;
z=quad(f,a,b,;
fprintf('z=%.4f\n',z)
z=
或:
f=inline('sqrt(3.*sin(t).*sin(t)+1)');
t=(0:
:
*pi;
y=f(t);
p=4*trapz(t,y);
fprintf('p=%.5f\n',p)
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- Matlab 软件 基础 数学 实验 程序 部分