常微分方程组边值.docx
《常微分方程组边值.docx》由会员分享，可在线阅读，更多相关《常微分方程组边值.docx（23页珍藏版）》请在冰豆网上搜索。
常微分方程组边值
常微分方程组边值问题解法
打靶法ShootingMethod（shooting.m）
%打靶法求常微分方程的边值问题
function[x,a,b,n]=shooting（fun,x0,xn,eps）
ifnargin<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'）
holdon
[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）
%边值常微分方程（组）函数
functionf=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）;
holdon
plot（x0,y0,’o’）
有限差分法FiniteDifferenceMethodsFDM（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））;
fori=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;
fori=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'）
holdon
[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,’-.’）
正交配置法OrthogonalCollocatioinMethodsCM
构造正交矩阵函数（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）
fori=1:
m
xm（i）=x0（m+1-i）;
end
xm（m+1）=1;
forj=1:
m+1
fori=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;
fori=2:
n+1
xn（i）=x1（n+2-i）;
end
xn（n+2）=1;
forj=1:
n+2
fori=1:
n+2
qn（j,i）=xn（j）^（i-1）;
ifj==0|i==1
cn（j,i）=0;
else
cn（j,i）=（i-1）*xn（j）^（i-2）;
end
ifj==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）;
%正交多项式求根（适用于对称问题）
functionp=symm（a,m,fm）%a为形状因子,m为配置点数,fm为权函数
fori=1:
m
c1=1;
c2=1;
c3=1;
c4=1;
forj=0:
i-1
c1=c1*（-m+j）;
iffm==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））;
%正交多项式求根（适用于非对称性问题）
functionp=unsymm（n,fn）
iffn==0
r
（1）=（-1）^n*n*（n+1）;%权函数为1
else
r
（1）=（-1）^n*n*（n+2）;%权函数为1-x
end
fori=1:
n-1
iffn==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
forj=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;
fori=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'）
holdon
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）
functionf=collfun1（x,y）
f
（1）=y
（2）;
f
（2）=4*y
（1）;
f=[f
（1）;f
（2）];
（collbc1.m）
functionf=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）
holdon
plot（sol.x,exp（2*sol.x）,'*'）真实解
精确解：
函数：
（collfun2.m）
functionf=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）
functionf=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）
holdon
plot（sol.x,（sol.x-1）.*exp（-sol.x）,'*'）真实解
精确解：
函数：
（collfun3.m）
functionf=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）
functionf=collbc3（a,b）
f=[2*a
（1）-a
（2）;b
（2）-1.5];
命令：
solinit=bvpinit（[1
升级会员 