How to solve a fourth order differential equation (boundary value problem) with constraint?

Hi everyone, I was trying to solve the following fourth order differential equation: A·x-B·x·w''''(x)=w(x) where A and B are constants,and 0<x<l, and the boundary conditions are w(0) = 0, w'(0) = 0, w''(l)=0, w'''(l)=0. I solved the equation with bvp4c. Now I want to add constraint to w''''(x), i.e., 0<w''''(x)<C, where C is a constant. How can I add this constraint? Thanks, JZ

Answer 1

[Torsten]

And as far as I see, a necessary and sufficient condition that w'''' is bounded on [0 l] is w(0) = 0, and that's what you already postulate as boundary condition for your problem.

[Jintao]

Sorry for the wrong boundary conditions, I've corrected it. So can I add this constraint? And how to do that?

Thanks a lot.

[Torsten]

As said, you don't need to add this constraint because you implicitly use it already: you postulate w(0) = 0 as boundary condition. In my opinion, this will ensure that w'''' is bounded:

Your equation reads

w'''' = (A*x-w)/(B*x)

Now this shows that x = 0 is the only critical point where w'''' could be unbounded. If you postulate that w(0) = 0, lim(x->0) w''''(x) = lim(x->0) (A*x-w)/(B*x) exists and equals (A-w'(0))/B by L'Hospital. Thus w'''' is bounded at 0 and thus on all of [0 l].

But I might be wrong: the time of mathematical studies is a long time ago.

What can be said is that you can't impose any additional constraints in bvp4c if they are not already implicit part of the boundary conditions. And if you think your differential equation could have a different solution with a constraint than it has without a constraint, a numerical solver like bvp4c is also not the suitable tool for investigation. Such a question needs theoretical analysis of the equation.

[Jintao]

Thanks

[Torsten]

Ok, I see now that you want to have w'''' bounded and positive. I'll think about the last point (which would mean that w(x) < A*x if B > 0 and w(x) > A*x if B < 0).

[Jintao]

What I am actually solving is the bending stress generated in an inclined steel rod (with one end encastrated) induced by wood swelling due to wetting. Based on Euler-Bernoulli beam theory, the fourth-order deflection of the steel rod (w''''(x)) is related to the distributed external load on the steel rod, and the magnitude of this external load should not exceed wood compressive yielding strength (C in the previous comments), that's why I want to add constraint to w''''(x). The following is my code,

clear all
L = 100; %STS penetration length
Xlow = 1e-6;
Xhigh = L;
MeshNumber = 100;
global ESTS %STS young's modulus (MPa)
ESTS = 226600;
global radius %STS radius
radius = 4;
global I %STS second moment of area
I = 0.25*pi*radius^4;
global alpha %swelling coefficient
alpha = 0.115*2;
global MCInitial %initial moisture content
MCInitial = 0.12;
global MCFinal %final moisture content
MCFinal = 0.30;
global theta %STS inclination angle
theta = 45;
global EWoodLongitudinal %Wood longitudinal young's modulus (MPa)
global EWoodRadial %Wood radial young's modulus (MPa)
if MCFinal == 0.12
    EWoodLongitudinal = 12400;
    EWoodRadial = 843.2;
elseif MCFinal == 0.16
    EWoodLongitudinal = 12400*(1-1.625*100*(MCFinal-MCInitial)/100);
    EWoodRadial = 843.2*(1-2.875*100*(MCFinal-MCInitial)/100);
elseif MCFinal == 0.21
    EWoodLongitudinal = 12400*(1-1.625*100*(MCFinal-MCInitial)/100);
    EWoodRadial = 843.2*(1-2.875*100*(MCFinal-MCInitial)/100);
elseif MCFinal == 0.30
    EWoodLongitudinal = 12400*(1-1.625*100*(MCFinal-MCInitial)/100);
    EWoodRadial = 843.2*(1-2.875*100*(MCFinal-MCInitial)/100);
else
    disp('Wrong final moisture content!!!Input agian!!!')
end
global EWoodTheta %Wood young's modulus from Hankinson fomula (MPa)
fenzi = EWoodLongitudinal*EWoodRadial;
fenmu = EWoodLongitudinal*sind(90-theta)^2+EWoodRadial*cosd(90-theta)^2;
EWoodTheta = fenzi/fenmu;
global beta %Wood compressive area reduction factor
beta = 0.5;
global b %Steel plate width
b = 1000;
options = bvpset('RelTol', 1e-6, 'AbsTol', 1e-6, 'NMax', 100000);
solinit = bvpinit(linspace(Xlow,Xhigh,MeshNumber),[0 0 0 0]);
sol = bvp4c(@mat4ode,@mat4bc,solinit,options);
x = linspace(Xlow,Xhigh,MeshNumber);
y = deval(sol,x);
plot(x,y(3,:)*ESTS*radius);

momentarray = y(3,:)*ESTS*I;
maxmoment = max(momentarray);
maxstress = maxmoment*radius/I;
minmoment = min(momentarray);
minstress = minmoment*radius/I;
function dxdy = mat4ode(x,y)
global ESTS
global I
global radius
global I
global alpha
global MCInitial
global MCFinal
global theta
global EWoodTheta
global beta
global b
%Coefficient A
A = alpha*(MCFinal-MCInitial)*sind(theta)*cosd(theta);
%Coefficient B
B = ESTS*I*tand(theta)/EWoodTheta/beta/(b+2*radius);
dxdy=[y(2)
    y(3)
    y(4)
    (A*x-y(1))/B/(0.5*(x*tand(theta)+sqrt(((b-2*radius)/2)^2+(x*tand(theta)+radius)^2)))]; %4th order equation
end
function res = mat4bc(ya,yb)
res=[ya(1) % w(0)=0 at x=0
    ya(2) % dw(x)/dx=0 at x=0
    yb(3) % M=0 at x=L
    yb(4)]; % Q=0 at x=L
end

[Torsten]

There is no way to impose an additional constraint apart from the 4 boundary conditions. If the solution of the differential equation you get does not satisfy your constraint and you think it should, then either the differential equation is wrong or implemented with errors or it has more than one solution.

[Jintao]

Thanks