微分方程:模型建立与求解微分方程
模型案例:最优化捕鱼问题
ode23,45
函数
dsolve
1.MATLAB求解函数
desolve('Dy = 2(1-0.01y)y','y(0) = 20','t')
result1=dsolve('Dy=2*(1-0.01*y)*y','y(0)=20','t') %符号型
syms y(t) %数值型
desolve(diff(y,t,1)==2*(1-0.01*y)*y,y(0)==20) %diff(y,t,1)--->y关于t的一阶导数
2.欧拉法与修正(不能求解解析解的情况下,数值解)
dy/dx = f(x,y) x>x0
区间取定离散点(步长一定)
取定步长: h,记 xn = x0 + nh, ( n = 1,2, ···, N )
称计算格式: yn+1 = yn + h f( xn, yn ) 为Euler公式。
常微分方程组的欧拉算法:
function testmain
close all
x0 = [0 0 0.1]';
tspan = 0:0.001:30;
[t,x] = ode45(@lorenz,tspan,x0);
subplot(1,2,1)
plot(t,x(:,1),'-',t,x(:,2),'*',t,x(:,3),'+')
subplot(1,2,2)
plot3(x(:,1),x(:,2),x(:,3)),grid on
set(gcf,'color','w')
figure
[T,Y]=myodesolver(@lorenz,t,x0);
plot3(Y(:,1),Y(:,2),Y(:,3),'r-')
function [T,Y]=myodesolver(odefun,T,y0)
% T时间变量离散化
% 要求y0为列向量,y0对应T(0)时刻的值
h = T(2)-T(1);%步长
Y = zeros(length(y0),length(T)); %最后转置
Y(:,1) = y0;
%y(1) = 10; % 对应t(1)时刻的值
for i= 2:length(T) % 使用欧拉法
Y(:,i) = Y(:,i-1) + h*odefun(T(i-1),Y(:,i-1));
end
Y = Y';
//
function testmain
close all
clc
t = 0:0.01:50;
h = 0.01;
b = 8/3;
c = 10;
r = 28;
x1(1) = 0;
x2(1) = 0;
x3(1) = 0.1;
for i= 2:length(t)
x1(i) = x1(i-1)+h*(-b*x1(i-1)+x2(i-1)*x3(i-1));
x2(i) = x2(i-1)+h*(-c*x2(i-1)+c*x3(i-1));
x3(i) = x3(i-1)+h*(-x1(i-1)*x2(i-1)+r*x2(i-1)-x3(i-1));
end
%y(i) = y(i-1) + h*dfun(t(i-1),y(i-1))
plot3(x1,x2,x3)
% [t,x] = ode45(@lorenz,[0,50],x0);
% subplot(2,1,1)
% plot(t,x(:,1),'-',t,x(:,2),'*',t,x(:,3),'+')
% subplot(2,1,2)
% plot3(x(:,1),x(:,2),x(:,3)),grid on set(gcf,'color','w')
function xdot = lorenz(t,x)
% 注意:xdot 必须是列向量
% x 也是列向量
b = 8/3; c = 10; r = 28;
xdot = [-b,0,x(2);0,-c,c;-x(2),r,-1]*x;
3.微分方程工具箱
ode23,ode45
ode(@f,自变量取值范围,自变量初值)
案例:核废料沉海是否符合安全规范:
案例叙述:
核废料沉海问题
——利用微分方程求沉底速度
摘 要
核电站处理废料时,常常将废料装入特制的圆桶中沉入海底,考虑到圆桶沉到海底时速度过大会撞破圆桶,则要计算出圆桶沉底的速度,若大于极限速度则圆桶会泄露,就不能在海中处理核废料;若小于则可以。
在假设的前提下,对圆桶进行合理的受力分析,列出牛顿第二定律的微分方程,解出结果v=45ft/s,速度v>40ft /s可知在海深300ft的地方投此圆桶,会在沉底时泄漏。
在整个解题时要注意:第一、我们只知道海深,并不知道圆桶沉底的时间,列微分方程时,v是关于圆桶里海平面变化的微分量 ;第二、单位的换算。
一、问题的提出
某核电站将核废料装于圆桶投入水深为300ft的海中,若圆桶在海中下降的速度大于或等于40ft/s时,会发生核泄漏,这是一种严重的事故,如果是这样的话就不能在海中处理核废料。问此核电站能否在指定海域处理核废料。
相关数据:
(1)1ft=0.3038m
(2)圆桶的容积为55gal (1gal=3.785L)
(3)装满废料时圆桶的质量为 527.436lbf (1lbf=1.425kg)
(4)海水浮力 B=470.32lbf
(5)海水阻力 D=Cv C=0.08
二、问题分析
圆桶在下落时受到向下的重力、向上的浮力和海水阻力,三力使圆桶受到向下的合力(如下图所示),使圆桶不断向下加速,直至沉底。由于不断加速,可知在沉底时速度最大,所以只要比较沉底速度和极限速度即可。
三、模型假设
1、在圆桶下沉过程中,圆筒内的核废料不会发生使圆桶质量变化的反应;
2、在圆桶下沉过程中,圆桶不会受到除重力、浮力、海水阻力以外的任何力;
3、在圆桶下沉过程中,海水在任何位置的属性都一样;
4、在圆桶下沉过程及沉底时,都不会出现加速度向上的情况;
四、符号约定
S —— 海水深度;
V —— 圆桶容积;
v —— 下沉速度;
vmax —— 圆桶的极限速度(大于此速度沉底,圆桶会撞破);
M —— 装核废物时圆桶的重量;
B —— 海水浮力;
D —— 海水阻力;
v(x(t)) —— t时刻圆桶与海平面的距离。
五、模型的建立及求解
根据圆桶在下沉过程中的受力情况,得:
mg-D-B=mdv/dt, (1)
因为在题中只给出海的深度,没有给出沉底所用的时间,所以对 dv/dt进行转换,dv/dt=(dv/dx)(dx/dt)=(dv/dx)v,得
mg-Cv-B=m(dv/dx)v (2)
又M=527.436lbf, 1lbf=3.785l, g=9.8m/s2, C=0.08, B=470.32lbf,
代入(2)中,解得 v = 45.1ft/s 。
六、结果分析
由于解得速度v大于极限速度vmax ,则在这个地方下投装满核废料的圆桶是会泄露的,所以不能在这里投放核废料。
在整个模型的建立过程中进行了一系列的假设,使得建立的模型理想化了,在实际中圆桶的下落过程是比较复杂的,海水的密度在不同的深度是不同的,则圆桶受到的浮力会应密度的变幻而变化;海水本身也有水流,会导致圆筒落入海底时,并不在投点的正下方,故在计算时,水深会发生变化;海水会对圆桶有一定的腐蚀作用,致使圆桶的强度变小;还有圆桶有个体差异,既使计算值小于极限速度,若差距不大,仍要做慎重考虑。总之,在实际处理时,还有很多问题要考虑。
编写MATLAB脚本如下:(采用欧拉算法与ode函数工具箱)
运算结果如图所示
显然当海水深度为300ft时桶触底速度远大于要求的极限速度,此方案安全性有待商榷
(欧拉算法在求解微分方程的数值解时需要带入一定相应的初始值,如本例中初始值为:深度为零时的桶的速度同样为零,然而带入微分方程时结果无意义,此时可采用初始值加上一个无限小作为新的带回有意义的初值,如本例中将v(1)=1e-6作为初值。)