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

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?

[Alexander]

Because you have endpoint constraints this is a Two Point Boundary Value Problem. To my knowledge there is not a way to solve TPBVP's through ODE, and in fact in general there is no guarantee that there is ANY solution to a problem like this.

[Alexander]

ODE is used for initial value problems, not boundary value problems. (just wanted to clarify that)

Answer 1

[bym]

you can try

bvp4c()

Answer 2

[Walter Roberson]

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.