常微分方程组(边值).doc

.
常微分方程组边值问题解法

打靶法Shooting Method（shooting.m）
%打靶法求常微分方程的边值问题
function [x,a,b,n]=shooting(fun,x0,xn,eps)
if nargin<3
eps=1e-3;
end
x1=x0+rand;
[a,b]=ode45(fun,[0,10],[0,x0]');
c0=b(length(b),1);
[a,b]=ode45(fun,[0,10],[0,x1]');
c1=b(length(b),1);
x2=x1-(c1-xn)*(x1-x0)/(c1-c0);
n=1;
while (norm(c1-xn)>=eps & norm(x2-x1)>=eps)
x0=x1;x1=x2;
[a,b]=ode45(fun,[0,10],[0,x0]');
c0=b(length(b),1);
[a,b]=ode45(fun,[0,10],[0,x1]');
c1=b(length(b),1)
x2=x1-(c1-xn)*(x1-x0)/(c1-c0);
n=n+1;
end
x=x2;

应用打靶法求解下列边值问题：

解：将其转化为常微分方程组的初值问题

命令：
x0=[0:0.1:10];
y0=32*((cos(5)-1)/sin(5)*sin(x0/2)-cos(x0/2)+1); 真实解
plot(x0,y0,'r')
hold on
[x,y]=ode45('odebvp',[0,10],[0,2]');
plot(x,y(:,1))
[x,y]=ode45('odebvp',[0,10],[0,5]');
plot(x,y(:,1))
[x,y]=ode45('odebvp',[0,10],[0,8]');
plot(x,y(:,1))
[x,y]=ode45('odebvp',[0,10],[0,10]');
plot(x,y(:,1))

函数：（odebvp.m）
%边值常微分方程(组)函数
function f=odebvp(x,y)
f(1)=y(2);
f(2)=8-y(1)/4;
f=[f(1);f(2)];

命令：
[t,x,y,n]=shooting('odebvp',10,0,1e-3)

计算结果：（eps=0.001）
t=11.9524
plot(x,y(:,1))
x0=[0:1:10];
y0=32*((cos(5)-1)/sin(5)*sin(x0/2)-cos(x0/2)+1);
hold on
plot(x0,y0,’o’)

有限差分法Finite Difference Methods FDM（difference.m）
同上例：

若划分为10个区间，则：

函数：（difference.m）
%有限差分法求常微分方程的边值问题
function [x,y]=difference(x0,xn,y0,yn,n)
h=(xn-x0)/n;
a=eye(n-1)*(-(2-h^2/4));
for i=1:n-2
a(i,i+1)=1;
a(i+1,i)=1;
end
b=ones(n-1,1)*8*h^2;
b(1)=b(1)-0;
b(n-1)=b(n-1)-0;
yy=a\b;
x(1)=x0;y(1)=y0;
for i=2:n
x(i)=x0+(i-1)*h;
y(i)=yy(i-1);
end
x(n)=xn;y(n)=yn;

命令：
[x,y]=difference(0,10,0,0,100);

计算结果：
x0=[0:0.1:10];
y0=32*((cos(5)-1)/sin(5)*sin(x0/2)-cos(x0/2)+1); 真实解
plot(x0,y0,'r')
hold on
[x,y]=difference(0,10,0,0,5);
plot(x,y,’.’)
[x,y]=difference(0,10,0,0,10);
plot(x,y,’--’)
[x,y]=difference(0,10,0,0,50);
plot(x,y,’-.’)

正交配置法Orthogonal Collocatioin Methods CM
构造正交矩阵函数（collmatrix.m）
%正交配置矩阵(均用矩阵法求对称性与非对称性正交配置矩阵)
function [am,bm,wm,an,bn,wn]=collmatrix(a,m,fm,n,fn)
x0=symm(a,m,fm); %a为形状因子;m为零点数;fm为对称的权函数(0为权函数1,非0为权函数1-x^2)
for i=1:m
xm(i)=x0(m+1-i);
end
xm(m+1)=1;
for j=1:m+1
for i=1:m+1
qm(j,i)=xm(j)^(2*i-2);
cm(j,i)=(2*i-2)*xm(j)^(2*i-3);
dm(j,i)=(2*i-2)*(2*i-3+(a-1))*xm(j)^(2*i-3+(a-1)-1-(a-1));
end
fmm(j)=1/(2*j-2+a);
end
am=cm*inv(qm);
bm=dm*inv(qm);
wm=fmm*inv(qm);
x1=unsymm(n,fn); %n为零点数;fn为非对称的权函数(0为权函数1,非0为权函数1-x)
xn(1)=0;
for i=2:n+1
xn(i)=x1(n+2-i);
end
xn(n+2)=1;
for j=1:n+2
for i=1:n+2
qn(j,i)=xn(j)^(i-1);
if j==0 | i==1
cn(j,i)=0;
else
cn(j,i)=(i-1)*xn(j)^(i-2);
end
if j==0 | i==1 | i==2
dn(j,i)=0;
else
dn(j,i)=(i-2)*(i-1)*xn(j)^(i-3);
end
end
fnn(j)=1/j;
end
an=cn*inv(qn);
bn=dn*inv(qn);
wn=fnn*inv(qn);

%正交多项式求根(适用于对称问题)
function p=symm(a,m,fm) %a为形状因子,m为配置点数,fm为权函数
for i=1:m
c1=1;
c2=1;
c3=1;
c4=1;
for j=0:i-1
c1=c1*(-m+j);
if fm==0
c2=c2*(m+a/2+j);%权函数为1
else
c2=c2*(m+a/2+j+1);%权函数为1-x^2
end
c3=c3*(a/2+j);
c4=c4*(1+j);
end
p(m+1-i)=c1*c2/c4/c3;
end
p(m+1)=1;%为多项式系数向量,求出根后对对称问题还应开方才是零点
p=sqrt(roots(p));

%正交多项式求根(适用于非对称性问题)
function p=unsymm(n,fn)
if fn==0
r(1)=(-1)^n*n*(n+1);%权函数为1
else
r(1)=(-1)^n*n*(n+2);%权函数为1-x
end
for i=1:n-1
if fn==0
r(i+1)=(n-i)*(i+n+1)*r(i)/(i+1)/(i+1);%权函数为1
else
r(i+1)=(n-i)*(i+n+2)*r(i)/(i+1)/(i+1);%权函数为1-x
end
end
for j=1:n
p(n+1-j)=(-1)^(j+1)*r(j);
end
p(n+1)=(-1)^(n+1);
p=roots(p);

应用正交配置法求解以下等温球形催化剂颗粒内反应物浓度分布，其浓度分布的数学模型为：

解：
（1）标准化
令，代入微分方程及边界条件得：

（2）离散化

（3）转化为代数方程组（以为例）

因为，所以整理上式得：

本例中的代数方程组为线性方程组，可采用线性方程组的求解方法；若为非线性方程组则采用相应的方法求解。

命令：N=3,权函数为1-x2
[am,bm,wm,an,bn,wn]=collmatrix(3,3,1,3,1);（只用对称性配置矩阵）
b1=bm;
for i=1:4
b1(i,i)=bm(i,i)-36;
end
a0=b1(1:4,1:3);
b0=-b1(1:4,4);
y=a0\b0;
y(4)=1;
p=exam31(3,3);（注意要对文件修改权函数为1-x2）
x=[0.3631,0.6772,0.8998,1]; %零点
plot(x,y,'o')
hold on
x0=0:0.1:1; %真实解
y0=sinh(6*x0)./x0/sinh(6);
plot(x0,y0,'r')

若权函数改为1，则以下语句修改，其他不变
[am,bm,wm,an,bn,wn]=collmatrix(3,3,0,3,1);（只用对称性配置矩阵）

p=exam31(3,3);（注意要对文件修改权函数为1）
x=[0.4058,0.7415,0.9491,1]; %零点

计算结果：
权函数为1- x2

权函数为1

边值问题的MatLab解法

精确解：

函数：（collfun1.m）
function f=collfun1(x,y)
f(1)=y(2);
f(2)=4*y(1);
f=[f(1);f(2)];
（collbc1.m）
function f=collbc1(a,b)
f=[a(1)-1;b(1)-exp(2)];

命令：
solinit=bvpinit([0:0.1:1],[1,1])
sol=bvp4c(@collfun1,@collbc1,solinit)
plot(sol.x,sol.y)
hold on
plot(sol.x,exp(2*sol.x),'*') 真实解

精确解：

函数：（collfun2.m）
function f=collfun2(x,y)
f(1)=y(2);
f(2)=(1-x.^2).*exp(-x)+2*y(1)-(x+1).*y(2);
f=[f(1);f(2)];
（collbc2.m）
function f=collbc2(a,b)
f=[a(2)-2;b(2)-exp(-1)];

命令：
solinit=bvpinit([0:0.1:1],[1,1]);
sol=bvp4c(@collfun2,@collbc2,solinit);
plot(sol.x,sol.y)
hold on
plot(sol.x,(sol.x-1).*exp(-sol.x),'*') 真实解

精确解：

函数：（collfun3.m）
function f=collfun3(x,y)
f(1)=y(2);
f(2)=(2-log(x))./x+y(1)./x-y(2).^2;
f=[f(1);f(2)];
（collbc3.m）
function f=collbc3(a,b)
f=[2*a(1)-a(2);b(2)-1.5];

命令：
solinit=bvpinit([1:0.1:2],[1,1]);
sol=bvp4c(@collfun3,@collbc3,solinit);
plot(sol.x,sol.y)
hold on
plot(sol.x,sol.x+log(sol.x),'*') 真实解

气流
在260的基础面上，为促进传热在此表面上增加纯铝的圆柱形肋片，其直径为25mm，高为150mm；该柱表面受到16气流的冷却，气流与肋片表面的对流传热系数为15，肋端绝热；肋片的导热系数为236，假设肋片的导热热阻与肋片表面的对流传热热阻相比可以忽略；试求肋片中的温度分布，及单个肋片的散热量为多少？

解：根据以上条件可知：导热热阻与对流热阻相比可以忽略，则在肋片径向上没有温度分布，在轴向上存在温度分布，外界气流与肋片的对流传热则可转化为内热源；故该问题为导热系数为常数的一维稳定热传导，其导热微分方程为：

边界条件为：时，（肋根）；时，（肋端绝热）。
分析解：，；传热量：
这是两点边值的常微分方程求解问题，故可转化为如下形式：

函数：（leipianfun.m leipianbc.m）
%圆柱形肋片(常微分方程组)
function f=leipianfun(x,y)
f(1)=y(2);
f(2)=15*pi*0.025/236/pi/0.025^2*4*(y(1)-16);
f=[f(1);f(2)];

%圆柱形肋片(边界条件)
function f=leipianbc(a,b)
f=[a(1)-260;b(2)];

命令：
solinit=bvpinit([0:0.01:0.150],[1,1]);
sol=bvp4c(@leipianfun,@leipianbc,solinit); %sol.y中每行对应sol.x节点的因变量值即：第一行为y1，第二行为y2值，依此类推；故第一行为函数值，第二行对应的一阶导数。
plot(sol.x,sol.y(1,:))
%以下为分析解
m=sqrt(15*pi*0.025/236/(pi/4*0.025^2))
c1=(exp(m*(sol.x-0.15))+exp(-m*(sol.x-0.15)))/2;
c2=(exp(m*(0.15))+exp(-m*(0.15)))/2;
t=16+(260-16)*c1/c2;
hold on
plot(sol.x,t,'*')
%计算传热量
q=-236*pi*0.025^2/4*sol.y(2,1);

计算结果：
q=40.1052W

在直径为20mm的圆管外安装环形肋片，其表面温度为260，肋片导热系数为45，置于16、对流传热系数为150的气流中；试根据单个环肋的传热量大小确定适宜的肋片高度和肋片厚度；并给出肋高为0.01m，肋厚为0.0003m环肋的温度分布。
解：近似为：导热系数为常数的一维稳定热传导，其导热微分方程为：

边界条件为：时，（肋根）；时，（肋端绝热）。
这是两点边值的常微分方程求解问题，故可转化为如下形式：

函数：（huanleifun.m huanleibc.m）
%环形肋片(常微分方程组)
function f=huanleifun(x,y,h,nada,delta,t0,tf)
f(1)=y(2);
f(2)=2*h/nada/delta*(y(1)-tf)-y(2)./x;
f=[f(1);f(2)];

%环形肋片(边界条件)
function f=huanleibc(a,b,h,nada,delta,t0,tf)
f=[a(1)-t0;b(2)];

命令：（hq.m）
%环肋不同肋高对散热量的影响
function [q,sol]=hq
h=150;
nada=45;
delta=0.0003;
t0=260;
tf=16;
r=0.01;
H=[0.001,0.002,0.003,0.004,0.005,0.006,0.007,0.008,0.009,0.01,0.012,0.014,0.016,0.018,0.02,0.03];
x=r+H;
for i=1:length(x)
solinit=bvpinit([r:0.001:x(i)],[1,1]);
sol=bvp4c(@huanleifun,@huanleibc,solinit,[],h,nada,delta,t0,tf);
q(i)=-nada*2*pi*r*delta*sol.y(2,1);
end

H=[0.001,0.002,0.003,0.004,0.005,0.006,0.007,0.008,0.009,0.01,0.012,0.014,0.016,0.018,0.02,0.03];
[q,sol]=hq;
plot(H,q)

计算结果：（肋厚为0.3mm）
由图可知，适宜的肋高可取0.005~0.015m。

命令：（dq.m）
%环肋不同肋厚对散热量的影响
function [q,sol]=dq
h=150;
nada=45;
delta=[0.0001,0.0002,0.0003,0.0004,0.0005,0.0006,0.0008,0.0009,0.001,0.0015,0.002,0.0025,0.003,0.0035,0.004,0.005];
t0=260;
tf=16;
r=0.01;
x=r+0.01;
for i=1:length(delta)
solinit=bvpinit([r:0.001:x],[1,1]);
sol=bvp4c(@huanleifun,@huanleibc,solinit,[],h,nada,delta(i),t0,tf);
q(i)=-nada*2*pi*r*delta(i)*sol.y(2,1);
end

delta=[0.0001,0.0002,0.0003,0.0004,0.0005,0.0006,0.0008,0.0009,0.001,0.0015,0.002,0.0025,0.
003,0.0035,0.004,0.005];
[q,sol]=dq;
plot(delta,q)

计算结果：（肋高为10mm）
由图可知，适宜的肋厚可取0.0005~0.002m，而在同一基础面上肋厚越大，则肋片数目越少，虽然肋厚增加散热量增大，但肋片数目减少导致的散热量减少幅度更大，所以适宜的肋厚综合考虑应取0.0005左右。

命令：（肋高为10mm，肋厚为0.3mm）
solinit=bvpinit([0.01:0.001:0.02],[1,1]);
sol=bvp4c(@huanleifun,@huanleibc,solinit,[],150,45,0.0003,260,16);
plot(sol.x-0.01,sol.y(1,:))

计算片状催化剂中等温一级不可逆反应的有效因子（等温内部效率因子）。已知数学模型为：
，
B.C. ，
其中为浓度；为催化剂的厚度；为Thiele模数，；；对片状催化剂，等温效率因子为：

解：本题为二阶ODE-BVP方程，转化为一阶ODEs：

B.C. ，

函数（efficiencyfun.m）
function f=efficiencyfun(x,y)
f(1)=y(2);
f(2)=36*y(1);
f=[f(1);f(2)];

函数（efficiencybcfun.m）
function f=efficiencybcfun(a,b)
f=[a(2);b(1)-1];

函数（efficiencybcfun.m）
function f=efficiency(x,xi,yi)
f=spline(xi,yi,x);

命令：
solinit=bvpinit([0:0.1:1],[1,1]);
sol=bvp4c(@efficiencyfun,@efficiencybcfun,solinit);
plot(sol.x,sol.y(1,:),':',sol.x,sol.y(2,:),'--')

quad(@efficiency,0,1,[1e-10],[],sol.x,sol.y(1,:))

结果：

0.1668 效率因子

教育资料