第四讲 bvp4c函数应用

为常数
j
bvp4c调用格式
sol bvp4c(@FunctionName,@BCFunction, solinit, options, p1, p2,L )
FunctionName函数的接口
function dydx=FunctionName(x,y,p1,p2,...) 其中x是一标量,y为 yj列向量, p1,p2,...等为fj 的参数，参数值已知，输出dydx是列向量， 对应于fj列向量
0 1, , M ,V分别为无量纲斜率，力矩及切力
梁的微分控制方程转化为一阶常微分方程组
y1 y
y3
d2y
d 2
y2
dy
d
y4
d3y
d 3
dy1
d
y2
dy2
d
y3
dy3
d
y4
dy4 q
d
用bvp4c函数求解受均布载荷的悬臂梁
y q 1
1 图1
悬臂梁的bvp4c函数求解
function CantileverBeam solinit=bvpinit(linspace(0,1,10),[0.5,0.5,0.5,0.5]); beamsol=bvp4c(@BeamODE, @BeamCantileverBC,solinit,[]); eta=linspace(0,1,50); y=deval(beamsol,eta); plot(eta,y(1,:))
BCFunction函数的创建
function Res=BCFunction (ya,yb,p1,p2,...) ya是表示 yj(a)的列向量,yb是表示 yj(b)的列 向量，即使不要求边界条件，已知参数p1,p2 等也必须出现在接口定义语句中，输出Res为 列向量。
调用bvp4c函数所涉及到的另外两个函数
EIy q(x)
工程梁的无量纲方程
如果定截面梁横向位移为w(x)，长为L，弹性模量是E，横截面的惯性 矩是I，单位长度载荷是 p0q(x),则无量纲方程如下：
q() d 4 y d 4
d P0L3 / EI
dy
d
M2
V
Vd P0 L
d3y d 3
其中 x / L, y y( ) w / h0，h0 P0L4 / EI，P0是常数，表示q的最大值
function bvpExample k=2; solinit=bvpinit(linspace(0,1,5),[-0.05,0.1]) exmpsol=bvp4c(@OdeBvp,@OdeBC,solinit,[],k); x=linspace(0,1,50); y=deval(exmpsol,x); plot(x,y(1,:)) function dydx=OdeBvp(x,y,k) dydx=[y(2);x-k*y(1)]; function res=OdeBC(ya,yb,k) res=[ya(1);yb(1)];
变量solinit为一结构？，可由函数bvpinit得到： solinit=bvpinit(x,y)，x向量为初始网格点的估计值， 向量y为每个yj的估计值。向量x和y的长度互不相关。
bvp4c的输出sol是一结构，为特定数量点对应的解。 为使曲线变得更光滑，需要在中间插入一些点，使 用：sxint=deval(sol,xint),xint为点向量，函数deval 根据这些点向量求解。sol为函数bvp4c的输出。
y 1 ( wl2 x wlx2 wx3 ) EI 2 2 6
y 1 ( wl2 x2 wlx3 wx4 ) EI 4 6 24
工程梁的微分控制方程推导
q M
M+dM
FP
dx
FN
dx
FN+d FN
FQ
FQ+dFQ
dM dx
FQ ,
dFQ q( x), dx
dFN p( x) dx
力图和弯矩图。
Q wl
+
x
M _
x
wl 2 2
求解过程
Q(x) w(l x)
M
(
x)
1 2
w(l
x)
2
用积分法，可推出梁的斜率（转角）和挠度表达式为
y
1 EI
(
wl2 x 2
wlx2 2
wx3 ) 6
c1
y
1 EI
(
wl 2 x2 4
wlx3 6
wx4 24
)
c1x
c2
利用x 0处的边界条件，y 0和y 0得c1 0和c2 0 最后梁的y和y的理论解为
求解结果显示
0 -0.01 -0.02 -0.03 -0.04 -0.05 -0.06 -0.07 -0.08 -0.09
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
solinit和exmpsol结构
solinit = x: [0 0.25000000000000 0.50000000000000
-0.035 -0.04 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
思考：用ode45函数和fsolve函数能用来求解梁的问题吗？
对于梁而言，本质是一个边值问题，梁的两端各为两个边 界条件，ode45只能求解初值问题，但对于梁的任一端，可 以先假设另两个初值已知，这样用ode45函数求可以求解控 制微分方程，但结果中包含两个未知量，然后再利用另一 端的已知边界条件，这样便建立了一个非线性方程组，可 用fsolve函数求得两个未知量。然后再用ode45求解。便可 得解。
function dydx = BeamODE(x, y) dydx = [y(2);y(3);y(4);1];
function bc=BeamCantileverBC(y0,y1) bc = [y0(1);y0(2);y1(3);y1(4)];
悬臂梁的bvp4c函数求解结果
0.14 0.12
0.1 0.08 0.06 0.04 0.02
结论：因此用ode45和fsolve联立是可以求解梁的问题！
bvp4c函数应用举例
求解下面这个单变量二阶常微分方程
d2y dx2
ky
x
0
x
1
边界条件为：
y(0) 0
y(1) 0
将上式重写为一个一阶方程组：
dy1 dx
y2
dy2 dx
x ky1
边界条件为：
y1 (0) 0 y 1(1) 0
y1 y dy
y2 dt
相关函数的创建
下面创建函数和子函数，实现一阶常微分方程的子函数为OdeBvp，记录边界条件 的函数为OdeBC,初始估计值由bvpinit给出，在0和1之间选择5个点，假定y1有一 个解-0.05，y2有一个解0.1,设k=2,则程序为:
0.75000000000000 1] y: [2x5 double] exmpsol = x: [1x7 double] y: [2x7 double] yp: [2x7 double] solver: 'bvp4c'
0
边值问题举例——工程梁
y w b
x
h
l
图1
材料力学解法回顾
根据材料力学中关于轴力，剪力和弯矩符号的规定：轴 力，拉力为正；剪力，绕隔离体顺时针方向转动者为正； 弯矩，使梁的下侧纤维受拉者为正。画出该悬臂梁的剪
简支梁(一端施加一额外力矩Mr=0.8）的bvp4c求解结果
0.014 0.012
0.01 0.008 0.006 0.004 0.002
0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0.005 0
-0.005 -0.01
-0.015 -0.02
-0.025 -0.03
0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
课堂练习—完成简支梁(一端施加一额外力矩Mr=0.8）的bvp4c函数程序
function SimpleBeam Mr=0.8; solinit=bvpinit(linspace(0,1,10),[0.5,0.5,0.5,0.5]); beamsol=bvp4c(@BeamODE, @SimpleBeamBC,solinit,[],Mr); eta=linspace(0,1,50); y=deval(beamsol,eta); plot(eta,y(1,:)) function dydx = BeamODE(x, y,Mr) dydx = [y(2);y(3);y(4);1]; function bc=SimpleBeamBC(y0,y1,Mr) bc = [y0(1);y0(3);y1(1);y1(3)-Mr];
第四讲 bvp4c函数的应用
bvp4c针对两点边界值问题进行数值求解，实现对 以下n个一阶常微分方程的数值解：
dy j dx
f j (x, y1, y2,L
, yn , q)
j=1,2,3L n
x的取值范围为a≤x≤b区间内，初始条件
yj (a) a
j 及y j (b) b
j,
j
1, 2L
, n, a j和b