常微分方程组边值Word格式.docx
《常微分方程组边值Word格式.docx》由会员分享，可在线阅读，更多相关《常微分方程组边值Word格式.docx（23页珍藏版）》请在冰豆网上搜索。
应用打靶法求解下列边值问题:
y100
解：
将其转化为常微分方程组的初值问题
命令:
xO=[O:
O.1:
1O];
真实解
y0=32*（（cos（5）-1）/sin（5）*sin（x0/2）-cos（x0/2）+1）;
plot（xO,yO,'
r'
）
holdon[x,y]=ode45（'
odebvp'
[0,10],[0,2]'
plot（x,y（:
1））
[x,y]=ode45（'
[0,10],[0,5]'
plot（x,y（:
[0,10],[0,8]'
[0,10],[0,10]'
函数：
（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（'
10,0,1e-3）
计算结果：
（eps=0.001）
t=11.9524
x0=[0:
1:
10];
holdon
plot（xO,yO,'
o'
（difference.m）
有限差分法FiniteDiffereneeMethodsFDM
同上例:
Yii2Yiyii
Yii2
若划分为10个区间，则：
h-YiYi18h2
4
h2
2——
1
2
Y1
Y2
8h2
Yn2
Yn1
.2
h
%有限差分法求常微分方程的边值问题function[x,y]=difference（xO,xn,yO,yn,n）h=（xn-xO）/n;
a=eye（n-1）*（-（2-hA2/4））;
fori=1:
n-2
a（i,i+1）=1;
a（i+1,i）=1;
b=ones（n-1,1）*8*hA2;
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-i）;
x（n）=xn;
y（n）=yn;
[x,y]=difference（0,10,0,0,100）;
plot（xO,yO,'
[x,y]=difference（0,10,0,0,5）;
plot（x,y,'
.'
[x,y]=difference（0,10,0,0,10）;
--'
[x,y]=difference（0,10,0,0,50）;
-.'
正交配置法OrthogonalCollocatioinMethodsCM
构造正交矩阵函数（collmatrix.m）
%E交配置矩阵（均用矩阵法求对称性与非对称性正交配置矩阵）
function[am,bm,wm,an,bn,wn]=collmatrix（a,m,fm,n,fn）
xO=symm（a,m,fm）;
%a为形状因子;
m为零点数;
fm为对称的权函数（0为权函数1,非0为权
函数1-xA2）
m
xm（i）=x0（m+1-i）;
xm（m+1）=1;
forj=1:
m+1
qm（j,i）=xm（j）A（2*i-2）;
cm（j,i）=（2*i-2）*xm（j）A（2*i-3）;
dm（j,i）=（2*i-2）*（2*i-3+（a-1））*xm（j）A（2*i-3+（a-1）-1-（a-1））;
fmm（j）=1/（2*j-2+a）;
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;
n+1
xn（i）=x1（n+2-i）;
xn（n+2）=1;
n+2
qn（j,i）=xn（j）A（i-1）;
ifj==0|i==1
cn（j,i）=0;
else
cn（j,i）=（i-1）*xn（j）A（i-2）;
ifj==0|i==1|i==2
dn（j,i）=O;
dn（j,i）=（i-2）*（i-1）*xn（jF（i-3）;
fnn（j）=1/j;
an=cn*inv（qn）;
bn=dn*inv（qn）;
wn=fnn*inv（qn）;
%E交多项式求根（适用于对称问题）
%E交多项式求根（适用于非对称性问题）
应用正交配置法求解以下等温球形催化剂颗粒内反应物浓度分布，其浓度分布的数学模型为：
（1）标准化
令xr/R，yC/Cs代入微分方程及边界条件得：
丄2x2dy36y
x2dxdx
x0単0
dx
x1,y1
（2）离散化
N1
Bjiyi36yj0
i1
j1,2,
（3）转化为代数方程组（以N3为例）
r1234yyyy
因为yNi目41，所以整理上式得:
B1136
B12
B13
B21
B2236
B23
B31
B32
B3336
B41
B42
B43
yiy2y3
B14
B24
B34
B4436
本例中的代数方程组为线性方程组，可采用线性方程组的求解方法；
若为非线性方程组则采用相应的方法求解。
N=3,权函数为1-x
[am,bm,wm,an,bn,wn]=collmatrix（3,3,1,3,1）;
b1=bm;
b1（i,i）=bm（i,i）-36;
a0=b1（1:
4,1:
3）;
b0=-b1（1:
4,4）;
y=aO\bO;
y（4）=1;
p=exam31（3,3）;
（注意要对文件修改权函数为x=[0.3631,0.6772,0.8998,1];
%零点
o'
x0=0:
0.1:
1;
%真实解
y0=sinh（6*x0）./x0/sinh（6）;
（只用对称性配置矩阵）
1-x2）
若权函数改为1,则以下语句修改，其他不变
[am,bm,wm,an,bn,wn]=collmatrix（3,3,0,3,1）;
（注意要对文件修改权函数为1）
x=[0.4058,0.7415,0.9491,1];
计算结果:
权函数为1
%
0.9
正交配置法
0.8
0.7
・
0.6
-
y0.5
0.4
0.3
X：
0.2
0.1
ii
1■
0.5
0.70.80.9
x
边值问题的MatLab解法
精确解:
2x
（collfuni.m）functionf=collfuni（x,y）f（i）=y
（2）;
f
（2）=4*y
（1）;
f=[f
（1）;
f
（2）];
（collbci.m）functionf=collbc1（a,b）f=[a
（1）-1;
b
（1）-exp
（2）];
solinit=bvpinit（[0:
1],[1,1]）sol=bvp4c（@collfun1,@collbc1,solinit）plot（sol.x,sol.y）
plot（sol.x,exp（2*sol.x）,'
*'
yiy2真实解
yx1y2y1x2ex0x1
y02,y11/e
yiy2
y21x2ex2y1
%02,y11
x1y2
1/e
（collfun2.m）
functionf=collfun2（x,y）
f
（2）=（1-x.A2）.*exp（-x）+2*y
（1）-（x+1）.*y
（2）;
f=[f
（1）；
（collbc2.m）
functionf=collbc2（a,b）f=[a
（2）-2;
b
（2）-exp（-1）];
solinit=bvpinit（[0:
1],[1,1]）;
sol=bvp4c（@co
升级会员 