常微分方程组(边值).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=1n-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
构造正交矩阵函数（collmatrixm）
%正交配置矩阵(均用矩阵法求对称性与非对称性正交配置矩阵)
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=1n+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=b11: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 fefficiencyfun(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 效率因子

教育资料