# Ordinary Differential Equation with boundary condition in z(0) and in z(L)

3 views (last 30 days)
Luca Amerio on 17 Oct 2011
Hi,
I'm quite new to MATLAB so excuse me if my question is very simple.
I'm trying to solve a ODE like the following one:
Ay'''' + By'' = C
with x in [0,L]
with the following boundary conditions:
• y(0)=0
• y'(0)=0
• y''(L)=0
• Ay'''(L) + By'(L) = C
Can someone suggest me a code to do this?
I've tried to look at the ODE help, but i can't figure how to impose the boundary condition on x=L.
EDIT: giving a look around i found the bvp4c function, but i'm having very hard times trying to get it working. Anyone has an hint for me?
##### 2 CommentsShowHide 1 older comment
Alexander on 17 Oct 2011
ODE is used for initial value problems, not boundary value problems. (just wanted to clarify that)

bym on 17 Oct 2011
you can try
bvp4c()

Walter Roberson on 18 Oct 2011
I was not able to get Maple's dsolve() to handle it directly. dsolve() can multiple handle boundary conditions usually, but in this case it complains about not being able to distinguish y(x) from y(L).
Anyhow, you can do this iteratively. Using Maple notation (MuPAD should be pretty similar):
Y := rhs(dsolve({A*(diff(y(x), x, x, x, x))+B*(diff(y(x), x, x)) = C, y(0) = 0, (D(y))(0) = 0}));
This yields
-_C3*A^(1/2)*sin(B^(1/2)*x/A^(1/2))/B^(1/2)-_C4*cos(B^(1/2)*x/A^(1/2))+(1/2)*C*x^2/B+_C3*x+_C4
where _C3 and _C4 are the two constants of integration to be determined.
We solve for _C3 using the next boundary condition, y``(L)=0:
C3 := simplify(solve(eval(diff(Y, x, x), x = L), _C3));
Then we take that back in to Y:
Y2 := simplify(subs(_C3=C3,Y));
This yields Y2 of
-(-(1/2)*A^(1/2)*(-2*B^(5/2)*cos(B^(1/2)*x/A^(1/2))*_C4+C*x^2*B^(3/2)+2*B^(5/2)*_C4)*sin(B^(1/2)*L/A^(1/2))+_C4*(-B^(5/2)*A^(1/2)*sin(B^(1/2)*x/A^(1/2))+x*B^3)*cos(B^(1/2)*L/A^(1/2))+C*(-B^(1/2)*A^(3/2)*sin(B^(1/2)*x/A^(1/2))+x*B*A))/(A^(1/2)*B^(5/2)*sin(B^(1/2)*L/A^(1/2)))
and now we solve for the final boundary condition:
C4 := simplify(solve(A*eval(diff(Y2,x,x,x),x=L) + B*eval(diff(Y2,x),x=L) = C, _C4));
This yields a C4 of
-(-B^(1/2)*A^(1/2)*(L-1)*sin(B^(1/2)*L/A^(1/2))+A)*C/(B^2*cos(B^(1/2)*L/A^(1/2)))
and this can then be substituted in to Y2 in order to get the complete y(x):
simplify(subs(_C4=C4,Y2));
which is
y(x) = ((A*B*(L-1)*sin(B^(1/2)*x/A^(1/2))-A^(1/2)*x*B^(3/2)*(-(1/2)*x-1+L))*cos(B^(1/2)*L/A^(1/2))-(A*B*(L-1)*sin(B^(1/2)*L/A^(1/2))-A^(3/2)*B^(1/2))*(-1+cos(B^(1/2)*x/A^(1/2))))*C/(B^(5/2)*A^(1/2)*cos(B^(1/2)*L/A^(1/2)))
This can be simplified slightly more by combining trig and moving factors, to get
y(x) = ((A^(1/2)*B^(1/2)*(L-1)*sin(B^(1/2)*x/A^(1/2))-x*(-(1/2)*x-1+L)*B)*cos(B^(1/2)*L/A^(1/2))+(-B^(1/2)*A^(1/2)*(L-1)*sin(B^(1/2)*L/A^(1/2))+A)*(-1+cos(B^(1/2)*x/A^(1/2))))*C/(B^2*cos(B^(1/2)*L/A^(1/2)))
You can get it a little more compact if you can assume that A and B are non-negative as the sqrt(B)/sqrt(A) can be merged to sqrt(B/A) and so on.