常微分方程组边值问题解法
打靶法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;
应用打靶法求解下列边值问题：
()()???
????==-
=010004822y y y dx
y d 解：将其转化为常微分方程组的初值问题
()?????
?
?????==-==t dx dy
y y y dx
dy y dx dy x 0011221
048
命令：
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 ） 同上例：
4822y dx y d -=?4822
11i i i i y h y y y -=+--+ 2121
842h y y h y i i i =+???
?
??--+- 若划分为10个区间，则：
????????????--=???????????????????????????
????????????
????? ??--???? ??--???
?
??--???? ??----08880842114211421142222212212
2
22h h h h y y y y h h h h n n M M O O O
函数：（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);
应用正交配置法求解以下等温球形催化剂颗粒内反应物浓度分布，其浓度分布的数学模型为：
???
?
?
???
?=====??? ??S
C C r dr dC r R C dr dC r dr d r
,10,03612
22 解：
（1）标准化
令R r x /=，S C C y /=代入微分方程及边界条件得：
???
?
?
???
?=====???
??1,10,036122y x dx dy x y
dx dy x dx d x
（2）离散化
0361
1
=-∑+=j N i i ji
y y B
1,2,1+=N j Λ
（3）转化为代数方程组（以3=N 为例）
???
???????=????????????????????----000036363636
432144434241
343332312423222114131211y y y y B B B B B B B B B B B B B B B B 因为141==+y y N ，所以整理上式得：
?
??
??
?????+----=????????????
??????---3636363644342414321434241333231232221131211B B B B y y y B B B B B B B B B B B B 本例中的代数方程组为线性方程组，可采用线性方程组的求解方法；若为非线性方程组则
采用相应的方法求解。
命令：N=3,权函数为1-x
2
[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-x 2
） 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
y 正交配置法真实解
边值问题的MatLab 解法
()()???==<<=''21,10104e y y x y y ? ()()???
??==='='2111221
1,104e
y y y y y y 精确解：
x e y 2=
函数：（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),'*') 真实解
()()
()()??
?='='<<-=-'++''-e y y x e x y y x y x /11,20101212 ? ()
()()()??
?
??
='='+-+-='='-e y y y x y e x y y y x /11,2012111
212221 精确解：
()x e x y --=1
函数：（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),'*') 真实解
()()()()?????='='-<<-=-'+''2/32,011221ln 212y y y x x x
x y y y ? ()()()()???????
==--+-='='2
/32,0112ln 212212
21221y y y y x y x x y y y 精确解：
x x y ln +=
函数：（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),'*') 真实解
在260C ?的基础面上，为促进传热在此表面上增加纯铝的圆柱形肋片，其直径为25mm ，高为150mm ；该柱表面受到16C ?气流的冷却，气流与肋片表面的对流传热系数为15K m W ?2/，肋端绝热；
肋片的导热系数为236K m W ?/，假设肋片的导热热阻与肋片表面的对流传热热阻相比可以忽略；试求肋片中的温度分布，及单个肋片的散热量为多
少？
解：根据以上条件可知：导热热阻与对流热阻相比可以忽略，则在肋片径向上没有温度分布，在轴向上存在温度分布，外界气流与肋片的对流传热则可转化为内热源；故该问题为导热系数为常数的一维稳定热传导，其导热微分方程为：
()C A t t hp dx
t d λλ∞-=Φ=&2
2 边界条件为：0=x 时，C 2600?=t （肋根）；H x =时，
0==H
x dx
dt
（肋端绝热）。
分析解：()
()[]()
mH ch H x m ch t t t t --+=∞∞0，C A hp m λ=；传热量：0
=-=x C
dx dt
A Q λ
这是两点边值的常微分方程求解问题，故可转化为如下形式：
()()()??
??
??
?≤≤==-='='∞H
x H y y A y y hp y y y C 0,0,260021221
λ
260?
函数：（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 的圆管外安装环形肋片，其表面温度为260C ?，肋片导热系数为45K W/m ?，置于16C ?、对流传热系数为150K W/m 2?的气流中；试根据单个环肋的传热量大小确定适宜的肋片高度和肋片厚度；并给出肋高为0.01m ，肋厚为0.0003m 环肋的温度分布。
解：近似为：导热系数为常数的一维稳定热传导，其导热微分方程为：
()()λδ
λλ∞∞-=-=
Φ
=??? ??t t h A t t hp dr dt r dr d r C 21& 边界条件为：1r r =时，C 2600?=t （肋根）；H r r r +==12时，
02
==r r dr
dt
（肋端绝热）。
这是两点边值的常微分方程求解问题，故可转化为如下形式：
()()()??
??
??
?
≤≤==--='='∞2
122112221,0,2602r x r r y r y x y y y h y y y λδ
函数：（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.0 16,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.0 16,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,:))
计算片状催化剂中等温一级不可逆反应的有效因子（等温内部效率因子）。已知数学模型为：
)(2
2
2c f dx
c d φ=，10< 00 ==x dx dc ，()11=c 其中c 为浓度；x 为催化剂的厚度；φ为Thiele 模数，6=φ；()c c f =；对片状催化剂，等温效率因子为： ()()()()()?=??=??=1 0101 0101 01dx c f dx c dx c f dx c f dx c f S η 解：本题为二阶ODE-BVP 方程，转化为一阶ODEs ： ()???='='12 2 21 y f y y y φ B.C. ()002=y ，()111=y 目录 引言 (1) 1预备知识 (2) 定义1.1（奇异Sturm-Liouville边值问题的正解） (2) 引理1.1.1 (2) 定义1.2（凸集的概念) (3) 定义1.3锥的定义 (3) 定义1.4（全连续算子的概念） (3) 1.5 （常微分边值问题的定义） (4) 定义1.6混合单调算子得定义） (4) 2 常微分方程边值问题正解得存在性 (5) 2.1 奇异Sturm-Liouville常微分边值问题的正解存在在 (5) 子 (8) 2.2 一类二阶边值问题的存在性 (9) 3一类混合单调算子应用 (11) 3.1一类混合单调算子的存在唯一性?........................ 错误!未定义书签。 3.2 求常微分边值问题的例题 (13) 结束语 (15) 参考文献 (15) 致 (16) 常微分方程边值问题与不动点定 （数学与统计学院 11级数学与应用数学2班）指导教师：攀峰 引言 从历史上看在有了微积分这个概念以后，紧接着出现了常微分方程。发展初期是属于“求通解”得时代，当人们从初期的热潮中结束要从维尔证明了卡帝方程中是一定不会存在一般性的初等解的时候开始的，并且柯西紧接着又提出了初值问题，常微分方程开始从重视“求通解”转向重视“求定解”的历史时代。 大学我们都学习了常微分方程这门学科，如果要研究它的定解问题，我们首先就会知道是常微分方程的初值问题。然而，在科学技术、生产实际问题中，我们还是提出了另一类定解问题-边值问题。对于常微分方程边值问题，伟大的科学家最早在解决二阶线性微分方程时，提出了分离变量法。[]1.在牛顿时期，科学家们已经提出过常微分的边值问题，牛顿也对常微分边值问题进行过研究，并且在1666年10月牛顿已经在这个领域取得了很大的成就，但是由于种种原因当时并没有整理成论文，所以没有及时出版。但在1687年他终于把在常微分方程上研究的成果发表了，虽然不是在数学著作中，却是他的一本力学著作中（《自然哲学的数学原理》）。 在微积分刚创立时期，雅克.伯努利来自瑞士的科学家提出了远著文明的问题-悬链线问题，紧着的地二年著名数学家莱布尼兹就给出了正确的解答，通过对绳子上个点受力分析，建立了以下方程 这个方程满足的定解条件是y（a）=α；y(b)=β.这是一个典型的常微分方程的边值问题。从这开始，常微分边值问题已经是科学家研究微分方程是不可或缺的工具，我就简单列举几个例子：（比如种族的生态系统；梁的非线性震动）等。对于怎么研究它， 第8章 常微分方程边值问题的数值解法 引 言 第7章介绍了求解常微分方程初值问题的常用的数值方法；本章将介绍常微分方程的边值问题的数值方法。 只含边界条件(boundary-value condition)作为定解条件的常微分方程求解问题称为常微分方程的边值问题(boundary-value problem). 为简明起见，我们以二阶边值问题为 则边值问题(8.1.1)有唯一解。 推论 若线性边值问题 ()()()()()(),, (),()y x p x y x q x y x f x a x b y a y b αβ'''=++≤≤?? ==? (8.1.2) 满足 (1) (),()p x q x 和()f x 在[,]a b 上连续； (2) 在[,]a b 上, ()0q x >, 则边值问题(8.1.1)有唯一解。 求边值问题的近似解，有三类基本方法： (1) 差分法(difference method)，也就是用差商代替微分方程及边界条件中的导数，最终化为代数方程求解； (2) 有限元法(finite element method)； (3) 把边值问题转化为初值问题，然后用求初值问题的方法求解。 差分法 8.2.1 一类特殊类型二阶线性常微分方程的边值问题的差分法 设二阶线性常微分方程的边值问题为 (8.2.1)(8.2.2) ()()()(),,(),(), y x q x y x f x a x b y a y b αβ''-=< 常微分方程组边值问题解法 打靶法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; 应用打靶法求解下列边值问题： ()()??? ????==- =010004822y y y dx y d 解：将其转化为常微分方程组的初值问题 ()????? ? ?????==-==t dx dy y y y dx dy y dx dy x 0011221 048 命令： 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)) 第十一章 常微分方程边值问题的数值解法 工程技术与科学实验中提出的大量问题是常微分方程边值问题．本章将研究常微分方程边值问题的数值求解方法.主要介绍三种边界条件下的定解问题和两大类求解边值问题的数值方法,打靶法算法和有限差分方法. 11.1 引言 在很多实际问题中都会遇到求解常微分方程边值问题. 考虑如下形式的二阶常微分方程 ),,(y y x f y '='', b x a <<, (11.1.1) 在如下三种边界条件下的定解问题: 第一种边界条件: α=)(a y , β=)(b y (11.1.2) 第二种边界条件: α=')(a y , β=')(b y (11.1.2) 第三种边界条件: ? ? ?=-'=-'101 0)()()()(b b y b y a a y a y βα, (11.1.13) 其中0 0, ,00000>+≥≥b a b a . 常微分方程边值问题有很多不同解法, 本书仅介绍打靶方法和有限差分方法. 11.2 打靶法 对于二阶非线性边值问题 ()()().,,βα==≤≤'=''b y a y b x a y y x f y ，，， (11.2.1) 打靶法近似于使用初值求解的情况． 我们需要利用一个如下形式问题初值解的序列： ()()v a w a w b x a w w x f w ='=≤≤'='')(,,，，，α， （11.2.2） 引进参数v 以近似原边界值问题的解．选择参数k v v =，以使： ()()β==∞ →b y v b w k k ,lim ， （11.2.3） 其中),(k v x w 定义为初值问题（11.2.2）在k v v =时的解，同时()x y 定义为边值问题（11.2.1）的解． 首先定义参数0v ，沿着如下初值问题解的曲线，可以求出点),(αa 对应的初始正视图 ()()v a w a w b x a w w x f w ='=≤≤'='')(,,，，，α． （11.2.4） 如果),(0v b w 不严格收敛于β，那么我们选择1v 等值以修正近似值，直到),(0v b w 严格逼近β． 为了取得合适的参数k v ，现在假定边值问题（11.2.1）有唯一解，如果),(v x w 定义为初始问题（11.2.2）的解，那么v 可由下式确定： 0),(=-βv b w ． （11.2.5） 由于这是一个非线性方程，我们可以利用Newton 法求解．首先选择初始值0v ，然后由下式生成序列 ),)(()),((111----- =k k k k v b dv dw v b w v v β，此处),(),)(( 11--=k k v b dv dw v b dv dw ， （11.2.6） 同时要求求得),)(( 1-k v b dv dw ，因为),(v b w 的表达式未知，所以求解这个有一点难度；我们只能得到这么一系列的值。 ，，，),(),(),(),(1210-??k v b w v b w v b w v b w 假如我们如下改写初值问题（11.2.2），使其强调解对x 和v 的依赖性 ()()v v a w v a w b x a v x w v x w x f w ='=≤≤'=''),(,),(),,(,，，，α，（11.2.7） 保留初始记号以显式与x 的微分相关．既然要求当k v v =时),)((v b dv dw 的值，那么我们需要求出表达式（11.2.7）关于v 的偏导数．过程如下： )),(),,(,(),(v x w v x w x v f v x v w '??=?''? ),()),(),,(,()),(),,(,(v x v w v x w v x w x w f v x v x w v x w x x f ??'??+??'??= ) ,()),(),,(,(v x v w v x w v x w x w f ?'?''??+ 又因为x 跟v 相互独立，所以当b x a ≤≤上式如下； 常微分方程组边值问题解法 打靶法Shooti ng Method (shoot in g.m ) % 丁靶法求常微分方程的边值问题 function [x,a,b ,n]=shooti ng(fu n, xO,x n, eps) if nargin<3 eps=1e-3; end x1=x0+ra nd; [a,b]=ode45(fu n, [0,10],[0,x0]'); c0=b(le ngth(b),1); [a,b]=ode45(fu n, [0,10],[0,x1]'); c1=b(le ngth(b),1); x2=x1-(c1-x n)*(x1-x0)/(c1-c0); n=1; while (no rm(c1-x n)>=eps & no rm(x2-x1)>=eps) x0=x1;x 仁x2; [a,b]=ode45(fu n,[ 0,10],[0,x0]'); cO=b(le ngth(b),1); [a,b]=ode45(fu n,[ 0,10],[0,x1]'); c1= b(le ngth(b),1) x2=x1-(c1-x n)*(x1-x0)/(c1-c0); n=n+1; end x=x2; 应用打靶法求解下列边值问题: y 10 0 解：将其转化为常微分方程组的初值问题 命令: xO=[O:O.1:1O]; y0=32*((cos(5)-1)/si n( 5)*si n(x0/2)-cos(x0/2)+1); plot(xO,yO,'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)) dy i dx y 2 dy 2 dx y i 0 y 4 y o dy dx X0 真实解 30 ' 12^4567^9 10 下载的，感觉不错，共享一下 常微分方程和常微分方程组的求解 一、实验目的： 熟悉Matlab 软件中关于求解常微分方程和常微分方程组的各种命令，掌握利用Matlab 软件进行常微分方程和常微分方程组的求解。 二、相关知识 在MATLAB 中，由函数dsolve()解决常微分方程（组）的求解问题，其具体格式如下： X=dsolve(‘eqn1’,’eqn2’,…) 函数dsolve 用来解符号常微分方程、方程组，如果没有初始条件，则求出通解，如果有初始条件，则求出特解。 例1：求解常微分方程1dy dx x y = +的MATLAB 程序为：dsolve('Dy=1/(x+y)','x')， 注意，系统缺省的自变量为t ，因此这里要把自变量写明。 结果为：-lambertw(-C1*exp(-x-1))-x-1 其中：Y=lambertw(X)表示函数关系Y*exp(Y)=X 。 例2：求解常微分方程2 '''0yy y -=的MATLAB 程序为： Y2=dsolve('y*D2y-Dy^2=0’,’x’) 结果为： Y2 =[ exp((x+C2)/C1)] [ C2] 我们看到有两个解，其中一个是常数。 例3：求常微分方程组253t t dx x y e dt dy x y e dt ?++=??? ?--=??通解的MATLAB 程序为： [X,Y]=dsolve('Dx+5*x+y=exp(t),Dy-x-3*y=exp(2*t)','t') 例4：求常微分方程组020 210cos ,224,0 t t t dx dy x t x dt dt dx dy y e y dt dt =-=?+-==??? ?++==??通解的MATLAB 程序为： [X,Y]=dsolve('Dx+2*x-Dy=10*cos(t),Dx+Dy+2*y=4*exp(-2*t)','x(0)=2','y(0)=0') 微分方程求解是系统仿真、数学模型实现以及很多工程问题求解的核心部分，应用MATLAB可以方便地对一阶常微分方程组进行求解，这里将对其基本方法进行介绍。值得注意的是，高阶微分方程组可以通过引进参变量化为一阶常微分方程组，也可以同样方便解决。 若有一个微分方程（组）的参变量为列向量,即,且它参变量随时间变化的微分方程可以有以下方程描述： 这里的f函数是一个列向量，即, i=1,2,3…,n，它可以是任意非线性函数。 则一般微分方程可以如此求解： [t,x]=ode45(f,timespan,x0) 对于刚性方程，即一些解变化缓慢，一些解变化剧烈，且两者相差较为悬殊的这种方程，通常调用ode15s而非o de45进行求解。 例1： 解：编写function或者用匿名函数 表达f=y-2*x/y即可； function dy=f(t,y) dy=y-2*t/y; end 命令： t=[0,1];%y0=1; [x,y]=ode45('f',t,1);%注意 这里的x相当于自变量t plot(x,y,x,sqrt(1+2*x)),legend('数值解','解析解'); 可见求解效果不错。 例2、 解：编写function function dx=f(t,x)%返回值是列向量 dx=[-x(2)-x(3); x(1)+0.2*x(2); 0.2+(x(1)-5.7)*x(3)]; end 命令： t=[0,100]; y0=[0 0 0]';%注意是列向量 [x,y]=ode45('f',t,y0); plot(x,y)； 例3、 这是一个二阶微分方程组，可以引进变量,由此ODE可以化成如下形式 可以采用和例2相同的方法求解： function dx=f(t,x) dx=[x(2); -(x(1)^2-1)*x(2)-x(1)]; End 二阶常微分方程边值问题的数值解法 摘要 求解微分方程数值解的方法是多种多样的，它本身已形成一个独立的研究方向，其要点是对微分方程定解问题进行离散化．本文以研究二阶常微分方程边值问题的数值解法为目标，综合所学相关知识和二阶常微分方程的相关理论，通过对此类方程的数值解法的研究，系统的复习并进一步加深对二阶常微分方成的数值解法的理解，为下一步更加深入的学习和研究奠定基础． 对于二阶常微分方程的边值问题，我们总结了两种常用的数值方法：打靶法和有限差分法．在本文中我们主要探讨关于有限差分法的数值解法．构造差分格式主要有两种途径：基于数值积分的构造方法和基于Taylor展开的构造方法．后一种更为灵活，它在构造差分格式的同时还可以得到关于截断误差的估计．在本文中对差分方法列出了详细的计算步骤和Matlab 程序代码，通过具体的算例对这种方法的优缺点进行了细致的比较．在第一章中，本文将系统地介绍二阶常微分方程和差分法的一些背景材料．在第二章中，本文将通过Taylor展开分别求得二阶常微分方程边值问题数值解的差分格式．在第三章中，在第二章的基础上利用Matlab求解具体算例，并进行误差分析． 关键词：常微分方程,边值问题，差分法，Taylor展开，数值解 The Numerical Solutions of Second-Order Ordinary Differential Equations with the Boundary Value Problems ABSTRACT The numerical solutions for solving differential equations are various. It formed an independent research branch. The key point is the discretization of the definite solution problems of differential equations. The goal of this paper is the numerical methods for solving second-order ordinary differential equations with the boundary value problems. This paper introduces the mathematics knowledge with the theory of finite difference. Through solving the problems, reviewing what have been learned systematically and understanding the ideas and methods of the finite difference method in a deeper layer, we can establish a foundation for the future learning. Matlab 求解常微分方程边值问题的方法：bvp4c 函数 常微分方程的边值问题，即boundary value problems ，简称BVP 问题，是指表达形式为 (,)((),())0'=??=?y f x y g y a y b 或(,,)((),(),)0'=??=? y f x y p g y a y b p 的方程组（p 是未知参数），在MA TLAB 中使用积分器bvp4c 来求解。 [命令函数] bvp4c [调用格式] sol=bvp4c(odefun,bcfun,solinit,options,p1,p2,…) sol 为一结构体，sol.x 、sol.y 、sol.yp 分别是所选择的网格点及其对应的y(x)与y'(x)数值; bvp4c 为带边值条件常微分方程积分器的函数命令；odefun 为描述微分方程组的函数文件；bcfun 为计算边界条件g(f(a),f(b),p)=0的函数文件；solinit 为一结构体，solinit.x 与solinit.y 分别是初始网格的有序节点与初始估计值，边界值条件分别对应a=solinit.x(l)和b=solinit.x(end)； options 为bvpset 命令设定的可选函数，可采用系统默认值；p1, p2…为未知参数。 例 求常微分方程0''+=y y 在(0)2=y 与(4)2=-y 时的数值解。 [解题过程] 仍使用常用方法改变方程的形式： 令1=y y ，21'=y y ，则原方程等价于标准形式的方程组1221 ?'=??'=-??y y y y ； 将其写为函数文件twoode.m ； 同时写出边界条件函数对应文件twobc.m ； 分别使用结构solinit 和命令bvp4c 确定y-x 的关系； 作出y-x 的关系曲线图。 [算例代码] solinit =bvpinit(linspace(0,4,5),[1 0]); % linspace(0,4,5)为初始网格，[1,0]为初始估计值 sol=bvp4c(@twoode,@twobc,solinit); % twoode 与twobc 分别为微分方程与边界条件的函数，solinit 为结构 x=linspace(0,4); %确定x 范围 y=deval(sol,x); %确定y 范围 plot(x,y(1,:)); %画出y-x 的图形 %定义twoode 函数（下述代码另存为工作目录下的twoode.m 文件） function dydx= twoode(x,y) %微分方程函数的定义 dydx =[y(2) -abs(y(1))]; %定义twobc 函数（下述代码另存为工作目录下的twobc.m 文件） function res= twobc(ya,yb); %边界条件函数的定义 res=[ya(1);yb(1)+2]; 常微分方程和常微分方程组的求解 一、实验目的： 熟悉Matlab 软件中关于求解常微分方程和常微分方程组的各种命令，掌握利用Matlab 软件进行常微分方程和常微分方程组的求解。 二、相关知识 在MATLAB 中，由函数dsolve()解决常微分方程（组）的求解问题，其具体格式如下： X=dsolve(‘eqn1’,’eqn2’,…) 函数dsolve 用来解符号常微分方程、方程组，如果没有初始条件，则求出通解，如果有初始条件，则求出特解。 例1：求解常微分方程1dy dx x y = +的MATLAB 程序为：dsolve('Dy=1/(x+y)','x')， 注意，系统缺省的自变量为t ，因此这里要把自变量写明。 结果为：-lambertw(-C1*exp(-x-1))-x-1 其中：Y=lambertw(X)表示函数关系Y*exp(Y)=X 。 例2：求解常微分方程 2 '''0yy y -=的MATLAB 程序为：Y2=dsolve('y*D2y-Dy^2=0’,’x’) 结果为： Y2 =[ exp((x+C2)/C1)] [ C2] 我们看到有两个解，其中一个是常数。 例3：求常微分方程组253t t dx x y e dt dy x y e dt ?++=??? ?--=??通解的MATLAB 程序为： [X,Y]=dsolve('Dx+5*x+y=exp(t),Dy-x-3*y=exp(2*t)','t') 例4：求常微分方程组020 210cos ,224,0 t t t dx dy x t x dt dt dx dy y e y dt dt =-=?+-==??? ?++==??通解的MATLAB 程序 为： 常微分方程论文 题目：常微分方程边值问题的数值解法 组长：数学132文洲 组员：数学131王琦 数学132姚瑶 信息132郭斌 院（系）：理学院 指导教师：岳宗敏 时间：2015年6月9日 常微分方程边值问题的数值解法 摘要：作为一类定解问题，补充条件由以自变量取某些值时，未知函数及其导数的值而定，称其为边值条件。许多物理和数学问题都归结为边值问题。本文介绍边值问题的待定常数法和格林函数。 关键词：边值问题 待定常数法 格林函数 Abstract: as a kind of definite solution problems, the supplementary conditions by took the certain values in the independent variable, the value of the unknown function and its derivative, referred to as boundary value conditions. Many physical and mathematical problems boil down to boundary value problems. In this paper, the boundary value problem of the method of undetermined constants and green's function. Keywords: boundary value problem Method of undetermined constants Green's function 11.1 引言 在很多实际问题中都会遇到求解常微分方程边值问题. 考虑如下形式的二阶常微分方程 ),,(y y x f y '='', b x a <<, (11.1.1) 在如下三种边界条件下的定解问题: 第一种边界条件: α=)a (y , β=)(b y (11.1.2) 第二种边界条件: α=')(a y , β=')(b y (11.1.2) 第三种边界条件: ?? ?=-'=-'1 01 0)()()()(b b y b y a a y a y βα, (11.1.13) 其中0 0, ,00000>+≥≥b a b a . 一、常微分方程组(ODEs) 简介 (1) 1. 简谐振动 (1) 2. 电路Vander Pol 方程 (1) 3. 生物种群的Volterra-Lotka 方程 (2) 4. 蝴蝶效应Lorenz 方程 (2) 二、MATLAB 数值求解ODEs的方法 (3) 1. 多变量常微分方程组的求解 (4) 2. 高阶常微分方程如何表示？ (4) 3. 相图和极限环怎么绘制？ (4) 个人在学习自动控制原理、现代控制理论、非线性动力学等课程时，经常遇到求解常微分方程组的问题。很多人知道MATLAB 是简便易行的一个工具，但是不会调用它自带的ode 求解器，往往还在自己编写单步欧拉法的程序，不仅求解精度差，而且程序不规范，还浪费了大量时间。以下我就工程中常见的一些非线性系统，利用MATLAB 自带的求解器，说明一下如何求解ODE 方程组、以及如何绘制相轨迹和极限环的问题。供相关专业工科大学生参考和借鉴。 一、常微分方程组(ODEs) 简介 以下列出了一些较为著名的非线性动力学系统的数学表达式，大都是由常微分方程组表达的。这种形式在工程中应用非常广泛，如力学中的非线性振动、航天领域的弹道计算、控制工程中的非线性系统等，由于自然界的大多数现象都表现出非线性，因此对于该种动力系统的研究以及微分方程的求解也具有重大的意义。以下列出一些工程应用中常见的一些由ODE 方程组所描述的动力系统。 1. 简谐振动 该式是一个2 阶非线性常微分方程。 2. 电路Vander Pol 方程 Fig 1. VanderPol 系统时域响应 3. 生物种群的 Volterra-Lotka 方程 Fig 2. Volterra-Lotka 方程时域响应（左） Fig 3 非线性动力学方程的极限环（右） 左图的捕食者 -猎物随时间变化的曲线表现出强烈的非线性，而状态变量 x 、y 的 变化却呈现出一个规则的鹅卵石状。 4. 蝴蝶效应 Lorenz 方程 微分方程边值问题的数值方法 本部分内容只介绍二阶常微分方程两点边值问题的的打靶法和差分法。 二阶常微分方程为 (,,),y f x y y a x b '''=≤≤ (1.1) 当(,,)f x y y '关于,y y '为线性时，即(,,)()()()f x y y p x y q x y r x ''=++，此时(1.1)变成线性微分方程 ()()(),y p x y q x y r x a x b '''--=≤≤ (1.2) 对于方程(1.1)或(1.2)，其边界条件有以下3类： 第一类边界条件为 (),()y a y b αβ== (1.3) 当0α=或者0β=时称为齐次的，否则称为非齐次的。 第二类边界条件为 (),()y a y b αβ''== (1.4) 当0α=或者0β=时称为齐次的，否则称为非齐次的。 第三类边界条件为 0101()(),()()y a y a y b y b ααββ''-=+= (1.5) 其中00000,0,0αβαβ≥≥+>，当10α=或者10β=称为齐次的，否则称为非齐次的。微分方程(1.1)或者(1.2)附加上第一类，第二类，第三类边界条件，分别称为第一，第二，第三边值问题。 1 打靶法介绍 下面以非线性方程的第一类边值问题(1.1)、(1.3)为例讨论打靶法，其基本原理是将边值问题转化为相应的初值问题求解。 【原理】假定()y a t '=，这里t 为解()y x 在x a =处的斜率，于是初值问题为 (,,) ()()y f x y y y a y a t α '''=?? =??'=? (1.6) 令z y '=，上述二阶方程转化为一阶方程组 一、常微分方程组(ODEs)简介 (1) 1. 简谐振动 (1) 2. 电路Vancfer Pol 方程 (1) 3. 生物种群的Volterra-Lotka方程 (2) 4. 蝴蝶效应Lorenz 方程 (2) 二、MATLAB数值求解ODEs的方法 (3) 1. 多变量常微分方程组的求解 (4) 2. 高阶常微分方程如何表示？ (4) 3. 相图和极限环怎么绘制？ (4) 个人在学习自动控制原理、现代控制理论、非线性动力学等课程时，经常遇到求 解常微分方程组的问题。很多人知道MATLAB是简便易行的一个工具，但是不会调用它自带的ode求解器，往往还在自己编写单步欧拉法的程序，不仅求解精度差，而且程序不规范，还浪费了大量时间。以下我就工程中常见的一些非线性系统，利用MATLAB自带的求解器，说明一下如何求解ODE方程组、以及如何绘制相轨迹和极限环的问题。供相关专业工科大学生参考和借鉴。 一、常微分方程组(ODEs)简介 以下列出了一些较为著名的非线性动力学系统的数学表达式，大都是由常微分方 程组表达的。这种形式在工程中应用非常广泛，如力学中的非线性振动、航天领 域的弹道计算、控制工程中的非线性系统等，由于自然界的大多数现象都表现出非线性，因此对于该种动力系统的研究以及微分方程的求解也具有重大的意义。以下列出一些工程应用中常见的一些由ODE方程组所描述的动力系统。 1. 简谐振动 x + y sinx = 0 该式是一个2阶非线性常微分方程。 2. 电路Vander Pol方程 x + p(x z—4-1 = 0 0 2 4 6 8 10 12 16 18 20 t/s Fig 1. VanderPol 系统时域响应 3. 生物种群的Volterra-Lotka 方程 =x （A - By ） =y （Cx 一 D ） 左图的捕食者-猎物随时间变化的曲线表现出强烈的非线性，而状态变量 X 、y 的 变化却呈现出一个规则的鹅卵石状。 4. 蝴蝶效应Lorenz 方程 x — — 10x 4 10y y = 28x — y - - xz 8 z = xy-^z 斗 Fig 2. Volterra-Lotka 方程时域响应（左） Fig 3非线性动力学方程的极限环（右） 摘要 本文主要研究二阶常微分方程边值问题的数值解法。对线性边值问题，我们总结了两类常用的数值方法，即打靶法和有限差分方法，对每种方法都列出了详细的计算步骤和Matlab程序代码，通过具体的算例对这两类方法的优缺点进行了细致的比较。 关键字：常微分方程边值问题；打靶法；差分法； ABSTRACT This article mainly discusses the numerical methods for solving Second-Order boundary value problems for Ordinary Differential Equations. On the one hand, we review two types of commonly used numerical methods for linear boundary value problems, i.e. shooting method and finite difference method. For each method, we give both the exact calculating steps , we compare the advantages and disadvantages in detail of these two methods through a specific numerical example. Key words：Boundary-Value Problems for Ordinary Differential Equations；Shooting Method；Finite Difference Method；常微分方程边值问题与不动点定论文
常微分方程边值问题的数值解法
常微分方程组(边值)
第十一章 常微分方程边值问题的数值解法汇总
常微分方程组(边值)
matlab常微分方程和常微分方程组的求解
常微分方程组的MATLAB求解范例
(完整版)二阶常微分方程边值问题的数值解法毕业论文
Matlab求解常微分方程边值问题的方法
matlab常微分方程和常微分方程组的求解
常微分方程边值问题的数值解法
常微分方程组的MATLAB求解方法
微分方程的边值问题
常微分方程组的MATLAB求解方法
二阶常微分方程边值问题的数值解法
相关主题
文本预览