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

Question

Jintao 2024-6-4

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2125476-how-to-solve-a-fourth-order-differential-equation-boundary-value-problem-with-constraint

评论： Jintao 2024-6-5

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2024-6-4

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2125476-how-to-solve-a-fourth-order-differential-equation-boundary-value-problem-with-constraint#answer_1467601

编辑：Torsten 2024-6-4

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.

7 个评论
显示 5更早的评论隐藏 5更早的评论

Jintao 2024-6-4

编辑：Torsten 2024-6-5

在 MATLAB Online 中打开

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 2024-6-5

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 2024-6-5