微分方程的MATLAB实现与欧拉算法

微分方程：模型建立与求解微分方程
模型案例：最优化捕鱼问题
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作为初值。）