数学建模综合实验放射性废料圆桶.doc

一.实验目的：1 归纳和学习求解常微分方程(组)的基本原理和方法； 2 掌握解析、数值解法，并学会用图形观察解的形态和进行解的定性分析； 3 MATLAB软件关于微分方程求解的各种命令。二.实验内容：1. 已知微分方程组 满足初始条件（1） 求上述微分方程组初值问题的特解(解析解)，并画出解函数的图形（2） 分别用 ode23、ode45 求上述微分方程组初值问题的数值解(近似解)，求解区间为利用画图来比较两种求解器之间的差异2分别用Euler 折线法和四阶 Runge-Kutta 法求解微分方程初值问题的数值解(步长h取0.1)，求解范围为区间0,3 3美国原子能委员会以往处理浓缩的放射性废料的方法，一直是把它们装入密封的圆桶里，然后扔到水深为90多米的海底。生态学家和科学家们表示担心，怕圆桶下沉到海底时与海底碰撞而发生破裂，从而造成核污染。原子能委员会分辨说这是不可能的。为此工程师们进行了碰撞实验，发现当圆桶下沉速度超过12.2m/s与海底相撞时，圆桶就可能发生碰裂。这样为避免圆桶碰裂，需要计算一下圆桶沉到海底时速度是多少？这时已知圆桶重量为239.46kg，体积为0.2058m3， 海水密度为1035.71 kg /m3。如果圆桶速度小于12.2m/s，就说明这种方法是安全可靠的，否则就要禁止用这种方法来处理放射性废料。假设水的阻力与速度大小成正比例，其正比例常数k=0.6。（1） 建立解决上述问题的微分方程数学模型。（2） 用数值和解析两种方法求解微分方程，并回答是否要禁止用这种方法来处理放射性废料。三. 实验方案（程序设计说明）开启软件平台MATLAB，开启MATLAB编辑窗口； 1. 根据微分方程求解步骤编写M文件 2. 保存文件并运行； 3. 观察运行结果(数值或图形)； 4. 根据观察到的结果和体会写出实验报告。 四. 实验步骤或程序（经调试后正确的源程序）第1题c1-1.m：Syms x y t x,y=dsolve(Dx+x+y=0,Dy+x-y=0, . x(0)=1, y(0)=0, t) x=simple(x) y=simple(y)ezplot(x,y,0,2);axis autoc1-2.m：函数文件vederpol.m :function xprime=verderpol(t,x)xprime=-x(1)-x(2);x(2)-x(1);命令文件vdpl.m:y0=1;0;t,x=ode23(verderpol,0,2,y0);x1=x(:,1);x2=x(:,2);plot(x1,x2,o-)fprintf(用ode23求得数值解为：)x=x1y=x2hold ont,x=ode45(verderpol,0,2,y0);x1=x(:,1);x2=x(:,2);plot(x1,x2,*)fprintf(用ode45求得数值解为：)x=x1y=x2第2题c2-1.m欧拉折线法clearf=sym(y-exp(x)*cos(x);a=0; b=3;h=0.1;n=(b-a)/h+1; % n=(b-a)/h;x=0; y=1;szj=x,y;for i=1:n-1 % i=1:n y=y+h*subs(f,x,y,x,y); x=x+h; szj=szj;x,y;endszjplot(szj(:,1),szj(:,2),-o)c2-2.m龙格-库塔法clear;f=sym(y-exp(x)*cos(x);a=0; b=3; h=0.1;n=(b-a)/h;x=0; y=1; szj=x,y;for i=1:n l1=subs(f,x,y,x,y); l2=subs(f,x,y,x+h/2,y+l1*h/2); l3=subs(f,x,y,x+h/2,y+l2*h/2); l4=subs(f,x,y,x+h,y+l3*h); y=y+h*(l1+2*l2+2*l3+l4)/6; x=x+h; szj=szj;x,y;endszjplot(szj(:,1),szj(:,2), dg-)第3题圆桶重量m=39.46千克 体积V=.2058立方米;阻力系数k=.6 下沉速度v下沉时的阻力D=kv; 下沉时间t(秒);圆桶在t时刻下沉深度y(t)(米) 水的密度p=1000重力加速度g=9.8运动物体受到的力等于该物体的质量与其运动加速度的乘积，即F=ma。所以有微分方程组解析法：c3-1.msyms m V p g k tt vvy=dsolve(m*D2y-m*g+p*V*g+k*Dy,y(0)=0,Dy(0)=0)y=subs(y,m,V,p,g,k,239.46,0.2058,1000,9.8,0.6);y=vpa(y,5);%求位移函数v=dsolve(m*Dv-m*g+p*V*g+k*v,v(0)=0)v=subs(v,m,V,p,g,k,239.46,0.2058,1000,9.8,0.6); v=vpa(v,5);%求速度函数s=y-90;tt=solve(s);%求到达海底90 米处的时间vv=subs(v,tt);%求到达海底90米处的速度vt=vv(1)数值法：c3-2.mfun=inline(9.8-(1035.71*9.8*0.2058)/239.46-0.6*v/239.46,v,t);t,v=ode23(fun,0,12.99939781354041734444456875868,0)plot (t,v,o-)ma