《数学建模综合实验放射性废料圆桶》由会员分享，可在线阅读，更多相关《数学建模综合实验放射性废料圆桶（8页珍藏版）》请在金锄头文库上搜索。
1、一.实验目的：1 归纳和学习求解常微分方程(组)的基本原理和方法； 2 掌握解析、数值解法，并学会用图形观察解的形态和进行解的定性分析； 3 MATLAB软件关于微分方程求解的各种命令。二.实验内容：1. 已知微分方程组 满足初始条件（1） 求上述微分方程组初值问题的特解(解析解)，并画出解函数的图形（2） 分别用 ode23、ode45 求上述微分方程组初值问题的数值解(近似解)，求解区间为利用画图来比较两种求解器之间的差异2分别用Euler 折线法和四阶 Runge-Kutta 法求解微分方程初值问题的数值解(步长h取0.1)，求解范围为区间0,3 3美国原子能委员会以往处理浓缩的放射性废
2、料的方法，一直是把它们装入密封的圆桶里，然后扔到水深为90多米的海底。生态学家和科学家们表示担心，怕圆桶下沉到海底时与海底碰撞而发生破裂，从而造成核污染。原子能委员会分辨说这是不可能的。为此工程师们进行了碰撞实验，发现当圆桶下沉速度超过12.2m/s与海底相撞时，圆桶就可能发生碰裂。这样为避免圆桶碰裂，需要计算一下圆桶沉到海底时速度是多少？这时已知圆桶重量为239.46kg，体积为0.2058m3， 海水密度为1035.71 kg /m3。如果圆桶速度小于12.2m/s，就说明这种方法是安全可靠的，否则就要禁止用这种方法来处理放射性废料。假设水的阻力与速度大小成正比例，其正比例常数k=0.6。
3、（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)ezplo
4、t(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求得数值解为：
5、)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=su
6、bs(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运动物体受到的力等于该物体
7、的质量与其运动加速度的乘积，即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 米处的时间
8、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);plot (t,v,o-)max_v=max(v)if max(v)12.2 fprintf(v(max)12.2,该方法不合理。n)else fprintf(v(max)12.2,则该种方法合理。)end五程序运行结果第1题c1-1.m：x=1/2*exp(2(1/2)*t)+1/4*2(1/2)*exp(-2(1/2)*t)-1/4*2(1/2)*exp(2(1/2)*t)+1/2
9、*exp(-2(1/2)*t) y =1/4*2(1/2)*exp(-2(1/2)*t)-1/4*2(1/2)*exp(2(1/2)*t)x=1/4*(2+2(1/2)*exp(-2*2(1/2)*t)-2(1/2)+2*exp(-2*2(1/2)*t)*exp(2(1/2)*t)y =1/4*2(1/2)*(exp(-2(1/2)*t)-exp(2(1/2)*t)c1-2.m用ode23求得数值解为：x = 1.0000 0.9999 0.9995 0.9975 0.9877.y = 0 -0.0001 -0.0005 -0.0025 -0.0125 用ode45求得数值解为：x = 1.0
10、000 0.9999 0.9999 0.9998 0.9998 y = 0 -0.0001 -0.0001 -0.0002 -0.0002 第2题：c2-1.mszj = 0 1.0000 0.1000 1.0000 0.2000 0.9900 0.3000 0.9693 0.4000 0.9373 0.5000 0.8936 0.6000 0.8383 0.7000 0.7718 0.8000 0.6949 0.9000 0.6093 1.0000 0.5174 1.1000 0.4223 1.2000 0.3282 1.3000 0.2407 1.4000 0.1666 1.5000 0.
11、1144 1.6000 0.0941 1.7000 0.1180 1.8000 0.2003 1.9000 0.3578 2.0000 0.6097 2.1000 0.9782 2.2000 1.4883 2.3000 2.1682 2.4000 3.0496 2.5000 4.1674 2.6000 5.5602 2.7000 7.2699 2.8000 9.3421 2.9000 11.8258 3.0000 14.7730c2-2.mszj = 0 1.0000 0.1000 0.9948 0.2000 0.9787 0.3000 0.9509 0.4000 0.9109 0.5000
12、0.8583 0.6000 0.7933 0.7000 0.7165 0.8000 0.6290 0.9000 0.5329 1.0000 0.4309 1.1000 0.3268 1.2000 0.2256 1.3000 0.1337 1.4000 0.0590 1.5000 0.0112 1.6000 0.0021 1.7000 0.0456 1.8000 0.1582 1.9000 0.3590 2.0000 0.6702 2.1000 1.1171 2.2000 1.7283 2.3000 2.5364 2.4000 3.5774 2.5000 4.8916 2.6000 6.5231 2.7000 8.5204 2.8000 10.9359 2.9000 13.8260 3.0000 17.2510第3题解析法：c3-1.my =g*(m-V*p)/k*t-m*g*(m-V*p)/k2+m*g*(m-V*p)/k2*exp(-k/m*t)v =g/k*m-g/k*V*p-exp(-k/m*t)*g*(m-V*p)/k圆桶到达海底的速度已超过12.2m/ s ，所以这
