y’’+ycos（x）=0,y（0）=1,y’（0）=0;
2．用向前欧拉公式和改进的欧拉公式求方程y’=y-2x/y,y（0）=1（0≤x≤1,h=0.1）的数值解，要求编写程序，并比较两种方法的计算结果，说明了什么问题？
4.Apollo卫星的运动轨迹的绘制
二、实验过程
1、
（1）
源程序：
y=dsolve（'Dy=y+2*x','y（0）=1','x'）
x=linspace（0,1,100）;
y=3*exp（x）-2*x-2;
plot（x,y）
实验结果：
y=3*exp（x）-2*x-2
1、
（2）
源程序:
y=dsolve（'D2y+y*cos（x）=0','y（0）=1','Dy（0）=0','x'）
实验结果：
Warning:
Explicitsolutioncouldnotbefound.
>Indsolveat194
y=[emptysym]
该方程没有解析解。
2、向前欧拉：
源程序：
euler.m:
function[x,y]=euler（fun,x0,xfinal,y0,n）;
ifnargin<5,n=50;
end
h=（xfinal-x0）/n;
x
（1）=x0;y
（1）=y0;
fori=1:
n
x（i+1）=x（i）+h;
y（i+1）=y（i）+h*feval（fun,x（i）,y（i））;
end
doty.m:
functionf=doty（x,y）;
f=y-2*x/y;
在matlab中输入：
[x,y]=euler（'doty',0,1,1,10）；
y1=sqrt（1+2.*x）;
plot（x,abs（y-y1）,'*'）%做出数值解与精确解的误差图
实验结果：
x=Columns1through10
00.10000.20000.30000.40000.50000.60000.70000.80000.9000
Column11
1.0000
y=Columns1through10
1.00001.10001.19181.27741.35821.43511.50901.58031.64981.7178
Column11
1.7848
改进欧拉：
源程序：
eulerpro.m:
function[x,y]=eulerpro（fun,x0,xfinal,y0,n）;
ifnargin<5,n=50;
end
h=（xfinal-x0）/n;x
（1）=x0;y
（1）=y0;
fori=1:
n
x（i+1）=x（i）+h;
y1=y（i）+h*feval（fun,x（i）,y（i））;
y2=y（i）+h*feval（fun,x（i+1）,y1）;
y（i+1）=（y1+y2）/2;
end
在matlab中输入：
[x,y]=eulerpro（'doty',0,1,1,10）；
plot（x,abs（y-y1）,'*'）
实验结果：
x=Columns1through10
00.10000.20000.30000.40000.50000.60000.70000.80000.9000
Column11
1.0000
y=Columns1through10
1.00001.09591.18411.26621.34341.41641.48601.55251.61651.6782
Column11
1.7379
实验结果分析：
由两张数值解与精确解的图像对比可得，改进欧拉公式误差更小。
4、源程序：
appollo.m:
functiondx=appollo（t,x）
mu=1/82.45;
mustar=1-mu;
r1=sqrt（（x
（1）+mu）^2+x（3）^2）;
r2=sqrt（（x
（1）-mustar）^2+x（3）^2）;
dx=[x
（2）
2*x（4）+x
（1）-mustar*（x
（1）+mu）/r1^3-mu*（x
（1）-mustar）/r2^3
x（4）
-2*x
（2）+x（3）-mustar*x（3）/r1^3-mu*x（3）/r2^3];
weixing.m:
x0=[1.2;0;0;-1.04935751];%x0（i）对应与xi的初值
options=odeset（'reltol',1e-8）;%该命令的另一种写法是options=odeset;options.reltol=1e-8;
tic
[t,y]=ode45（@appollo,[0,20],x0,options）;%t是时间点，y的第i列对应xi的值，t和y的行数相同
toc
plot（y（:
1）,y（:
3））%绘制x1和x3，也就是x和y的图形
title（'Appollo卫星运动轨迹'）
xlabel（'X'）
ylabel（'Y'）
实验结果：
Elapsedtimeis0.085081seconds.
应用实验
一、问题重述
8.两种生物种群竞争模型
两种相似的群体之间为了争夺有限的同一种食物来源和生活空间而进行生存竞争时，往往是竞争力较弱的种群灭亡，而竞争力较强的种群达到环境容许的最大数量。
假设有甲、乙两个生物种群，当它们各自生存于一个自然环境中，均服从Logistic规律，即有
其中x1（t）,x2（t）分别为两种生物种群在时刻t的数量，λ1，λ2分别为其自然增长率，N1，N2是它们各自的最大容量。
当两个种群在同一个自然环境下生存时，乙消耗的同一自然资源对甲的增长产生了阻滞作用，设为
甲对乙的阻滞作用设为
由于生物种群的数量很大，可视为时间t的连续可微函数。
生物种群的相互竞争模型为
1）1的含义是，对于供养甲的资源而言，单位数量乙（相对N2）的消耗为单位数量甲（相对N1）消耗的1倍，对2可作相应的解释。
计算x1（t）,x2（t）,画出图形及相轨迹图。
解释其解变化过程。
2）改变λ1，λ2,N1,N2,而α1，α2不变，计算并分析结果；若α1=1.5，α2=0.7，再分析结果。
由此能得到什么结论。
二、数学模型的建立与求解（一般应包括模型、求解步骤或思路，程序放在后面的附录中）
由于参数输入不便，所以以s和r替代，r1，r2为它们的固有增长率，n1，n2为它们的最大容量。
s1的含义是，对于供养甲的资源来说，单位数量的乙（相对n2）的消耗为单位数量甲（相对n1）消耗的s1倍，对s2可以作相应解释。
（1）设r1=r2=1,n1=n2=100，s1=0.5,s2=2，初值x0=y0=10，输入程序，做出图像。
图1x1（t）,x2（t）图像
图2x（y）图像
可见最后数值稳定在x=100,y=0上，即物种甲达到最大值，物种乙灭绝。
物种乙开始一段时间数量稍稍有所增长，10年后就渐渐灭绝了，最后稳定状态就只剩下甲物种。
（2）改变r1,r2（r1=r2=0.2）
改变后，甲乙两物种最终结果仍然是甲达到数量极限而乙灭绝，但与原图相比，变化减缓了。
（3）改变n1,n2:
n1=300,n2=120:
改变n并不影响两物种的结局，这条曲线更为夸张了。
（4）改变s:
s1=1.5，s2=0.7:
结果正和s1=0.5,s2=2时相反，最后甲物种灭绝，乙物种存活并达到数量极限。
如果这时改变r1,r2,n1,n2,x0,y0这些参数，变化趋势和上面列举的相同（甲乙相反），这从方程的对称性上可以求证。
所以可以得出结论，由s1,s2的物理意义，当某个s1或者s2大于1时（另一个小于1），它将严重消耗其作用的物种的生存资源，最终的结果是致使此物种灭绝。
三、附录（程序等）
fun.m:
functiondx=fun（t,x,r1,r2,n1,n2,s1,s2）
dx=[r1*x
（1）*（1-x
（1）/n1-s1*x
（2）/n2）;r2*x
（2）*（1-s2*x
（1）/n1-x
（2）/n2）];
s51.m:
h=0.1;%所取时间点间隔
ts=[0:
h:
30];%时间区间
x0=[10,10];%初始条件
opt=odeset（'reltol',1e-6,'abstol',1e-9）;%相对误差1e-6，绝对误差1e-9
[t,x]=ode45（@fun,ts,x0,opt,1,1,100,100,0.5,2）;%使用5级4阶龙格—库塔公式计算
plot（t,x,'.-'）,grid%输出x1（t）,x2（t）的图形
gtext（'x1（t）'）,gtext（'x2（t）'）,pause
plot（x（:
1）,x（:
2）,'.-'）,grid,%作相轨线
总结与体会
通过本次实验，我掌握了利用matlab解微分方程的方法，了解了数值解与解析解的区别，学会了新知识欧拉公式，更通过应用实验，加深了对其的理解和解实际问题的能力。
教师签名
年月日
升级会员 