重庆大学 数学实验 2方程求解Word格式.docx
《重庆大学 数学实验 2方程求解Word格式.docx》由会员分享,可在线阅读,更多相关《重庆大学 数学实验 2方程求解Word格式.docx(17页珍藏版)》请在冰豆网上搜索。
成绩
√
实验目的
[1]复习求解方程及方程组的基本原理和方法;
[2]掌握迭代算法;
[3]熟悉MATLAB软件编程环境;
掌握MATLAB编程语句(特别是循环、条件、控制等语句);
[4]通过范例展现求解实际问题的初步建模过程;
通过该实验的学习,复习和归纳方程求解或方程组求解的各种数值解法(简单迭代法、二分法、牛顿法、割线法等),初步了解数学建模过程。
这对于学生深入理解数学概念,掌握数学的思维方法,熟悉处理大量的工程计算问题的方法具有十分重要的意义。
基础实验
一、实验内容
1.方程求解和方程组的各种数值解法练习
2.直接使用MATLAB命令对方程和方程组进行求解练习
3.针对实际问题,试建立数学模型,并求解。
1.用图形放大法求解方程xsin(x)=1.并观察该方程有多少个根。
程序如下:
x=-5:
0.11:
5;
y=x.*sin(x)-1;
plot(x,y);
grid
取x=-1到-1.3可锝
取x=-1.15到-1.1可锝
可得解为x=-1.115
2.将方程x5+5x3-2x+1=0改写成各种等价的形式进行迭代,观察迭代是否收敛,并给出解释。
①迭代函数为
,算法设计为:
x1=0;
x2=(x1^5+5*x1^3+1)/2;
whileabs(x1-x2)>
10^(-5)
x1=x2;
x2=(x1^5+5*x1^3+1)/2;
end
x1
输出结果为:
x1=Inf
因此x=(x)迭代不收敛,则不直接使用(x)迭代,用加速迭代函数
,
算法设计为:
x2=(-4*x1^5-10*x1^3+1)/(-5*x1^4-15*x1^2+2);
x1=-0.7685
②迭代函数为
x1=1;
x2=((2*x1-x1^5-1)/5)^(1/3);
x2=((2*x1-x1^5-1)/5)^(1/3);
x1=Inf-Infi
x2=((0.4*x1-0.2*x1^5-0.2)^(1/3)-1/15*(0.4*x1-0.2*x1^5-0.2)^(-2/3)*(2*x1-5*x1^5))/(1-(1/15*(0.4*x1-0.2*x1^5-0.2)^(-2/3)*(2-5*x1^4)));
x1=0.4004+0.2860i
③迭代函数为
x2=(2*x1-5*x1^3-1)^(1/5);
fork=1:
100
x2=(2*x1-5*x1^3-1)^(1/5);
x1=2.0162-0.8223i
若用加速迭代函数
x2=((2*x1-5*x1^3-1)^(1/5)-1/5*(2*x1-5*x1^3-1)^(-4/5)*(2*x1-15*x1^3))/(1-1/5*(2*x1-5*x1^3-1)^(-4/5)*(2-15*x1^2));
x2=((2*x1-5*x1^3-1)^(1/5)-1/5*(2*x1-5*x1^3-1)^(-4/5)*(2*x1-15*x1^3))/(1-1/5*(2*x1-5*x1^3-1)^(-4/5)*(2-15*x1^2));
x1=-0.1483+0.7585i
④迭代函数为
x2=0.2*(2/x1-1/x1^2-x1^3);
x2=0.2*(2/x1-1/x1^2-x1^3);
输出结果为
x1=NaN
x2=((2/x1-1/x1^2-x1^3)-x*(-2/x1^2+2/x1^3-3*x1^2))/(5-(-2/x1^2+2/x1^3-3*x1^2));
x2=((2/x1-1/x1^2-x1^3)-x*(-2/x1^2+2/x1^3-3*x1^2))/(5-(-2/x1^2+2/x1^3-3*x1^2));
x1=3.4802308631248458912724395623836
⑤迭代函数为
x2=2/x1^3-5/x1-1/x1^4;
x2=2/x1^3-5/x1-1/x1^4;
x1=1.8933
x2=((2/x1^3-5/x1-1/x1^4)-x*(-6/x^4+5/x^2+4/x^5))/(1-(-6/x^4+5/x^2+4/x^5));
x2=((2/x1^3-5/x1-1/x1^4)-x*(-6/x^4+5/x^2+4/x^5))/(1-(-6/x^4+5/x^2+4/x^5));
x1=1.7968059417612661783255756706113
3.求解下列方程组
(1)
①用solve()对方程组求解,算法设计为:
[x1,x2]=solve('
2*x1-x2-exp(-x1)'
'
-x1+2*x2-exp(-x2)'
)
x1=.56714329040978387299996866221036
x2=.56714329040978387299996866221036
②用fsolve()对方程组求解:
建立名为fun1.m的M文件,算法设计如下:
functionf=fun1(x)
f
(1)=2*x
(1)-x
(2)-exp(-x
(1));
f
(2)=-x
(1)+2*x
(2)-exp(-x
(2));
在函数体外部调用此函数:
>
y=fsolve('
fun1'
[1,1],1)
Optimizationterminated:
first-orderoptimalityislessthanoptions.TolFun.
y=0.56710.5671
(2)①用solve()对方程组求解,算法设计为:
[x1,x2,x3]=solve('
x1^2-5*x2^2+7*x3^2+12'
3*x1*x2+x1*x3-11*x1'
2*x2*x3+40*x1'
double(x1)
double(x2)
double(x3)
ans=1.0e+002*
0.0100
0
-0.0031
-3.8701+0.3270i
-3.8701-0.3270i
ans=
5.0000
1.5492
-1.5492
2.9579
-0.3123-50.8065i
-0.3123+50.8065i
1.0e+002*
-0.0400
0+0.0131i
0-0.0131i
0.0213
0.1194+1.5242i
0.1194-1.5242i
建立名为fun2.m的M文件,算法设计如下:
functionf=fun2(x)
f
(1)=x
(1)^2-5*x
(2)^2+7*x(3)^2+12;
f
(2)=3*x
(1)*x
(2)+x
(1)*x(3)-11*x
(1);
f(3)=2*x
(2)*x(3)+40*x
(1);
fun2'
[1,1,1],1)
y=
0.00001.54920.0000
4.编写用二分法求方程根的函数M文件。
思路:
共建立两个M文件,第一个M文件供输入函数用,第二个M文件调用第一个M的运算值来执行二分法运算过程。
算法设计(以解方程x2-3x+2=0为例):
functionf=fun0(x)
f=x^2-3*x+2;
functionf=bisection(x)
m=x
(1);
n=x
(2);
while(n-m)>
iffun0(m)==0
f=m;
break;
elseiffun0(n)==0
f=n;
elseiffun0((m+n)/2)==0
f=(m+n)/2;
elseiffun0(m)*fun0((m+n)/2)<
n=(m+n)/2;
else
m=(m+n)/2;
end
在函数体外部调用函数,输入:
x=[0.25,3];
root=bisection(x)
输出结果
root=1.0000
综合实验
炮弹发射角的问题
炮弹发射视为斜抛运动,已知初始速度为200m/s,问要击中水平距离360m、垂直距离160m的目标,当忽略空气阻力时,发射角应多大?
此时炮弹的运行轨迹如何?
试进行动态模拟。
进一步思考:
如果要考虑水平方向的阻力,且设阻力与(水平方向)速度成正比,系数为0.1(1/s),结果又如何?
二、问题分析
炮击目标确定后,如何调整发射角度使炮弹能准确地落在目标位置处爆炸,需要考虑一下三个方面的变化量:
1.发射速度大小不变的情况下,发射角度是可以随意调整的;
2.对于某固定的发射角,炮弹在运动过程中,纵向分速度受重力加速度的影响而随时间不断地变化;
3.若考虑水平方向的阻力,则阻力带来的水平加速度也会对水平分速度造成随时间变化的影响。
三、数学模型的建立与求解
求解步骤或思路:
1.首先建立一个函数bomb1,实现以下功能:
对于每个输入的确定的水平分速度,都能计算出炮弹运动到横坐标为360时,与点(360,160)的距离;
2.再建立一个函数bomb2,实现以下功能:
对于每个输入的确定的水平分速度,都能描绘出炮弹的等时间间隔的轨迹;
3.最后建立一个M文件bomb0,利用循环不断改变水平分速度,调用函数bomb1和bomb2,实现以下功能:
寻找出炮弹运动到横坐标为360时,与点(360,160)的距离小于某限定值的水平分速度,最后通过这一求出的水平分速度,求出炮弹发射角度及描绘炮弹运动轨迹。
4.细节问题:
为了使炮弹的轨迹更连贯,时间间隔应该尽可能地小;
时间间隔减小,会增加运算次数,增加运算时间,因此有必要先选取合理计算范围,缩短计算时间;
选取合理计算范围方法:
先输入几个间隔较大的水平分速度值,估计范围,再逐步逼近范围。
四、实验结果及分析
无空气阻力情况下炮弹的运动图像及炮弹发射夹角:
angle=26.4244
有空气阻力情况下炮弹的运动图像及炮弹发射夹角:
angle=24.6241
五、附录(程序等)
无空气阻力:
bomb10.m:
x
(1)=178.8;
delta=0.0001;
f
(1)=bomb11(x
(1));
fmin=1000;
i=0;
whilex<
180
i=i+1;
f(i)=bomb11(x(i));
iff(i)<
fmin
fmin=f(i);
xmin=x(i);
x(i+1)=x(i)+delta;
angle=acos(xmin/200)/pi*180
bomb12(xmin)
bomb11.m:
functionf=bomb11(x)
Vy=sqrt(200^2-x^2);
A=0+1i*0;
Point=360+1i*160;
t=0;
while(x*t<
=360)
t=t+0.01;
Vy=Vy-9.8*0.01;
A=A+(x+1i*Vy)*0.01;
f=abs(A-Point);
bomb12.m:
functionbomb12(x)
holdon;
plot(x*t,sqrt(200^2-x^2)*t-4.9*t^2,'
.'
);
holdoff;
grid;
有空气阻力:
bomb20.m:
x
(1)=181;
f
(1)=bomb21(x
(1));
183
f(i)=bomb21(x(i));
bomb22(xmin)
bomb21.m:
functionf=bomb21(x)
Vx=x;
S=0;
while(S<
Vx=Vx-0.1*Vx*0.01;
A=A+(Vx+1i*Vy)*0.01;
S=S+Vx*0.01;
bomb22.m:
functionbomb22(x)
plot(S,sqrt(200^2-x^2)*t-4.9*t^2,'
年月日