solving simultaneous equations of fourth order PDE in MATLAB

I'm trying to solve this system of ordinary differential equation with the boundary conditions : w(0)=w(L)=u(0)=u(L)=w'(0)=w'(L)=u'(0)=u'(L)=0.

I'm new to MATLAB and am unable to solve it using ODE toolbox.

Can someone with knowledge on the topic help me with this?

[darova]

How to express and ?

[Devansh Gupta]

In terms of x and L. w and u are deflections in z and x direction respectively. Px and Pz are constant loads im x and z direction resp.

[darova]

I mean that highest range of ω and u should be expressed from the system to solve ODE

[Devansh Gupta]

I'm not sure what woud that be. It's a problem of a composite beam bending under a Uniformly Distributed Load. What could be the Range in that case?

[darova]

Here is the example. Assume you have this system of equations

To solve it you need to express highest x and y

And code for this will be

f = @(u) [u(2)  % x'
    1-u(4)      % x'' = 1 - y'
    u(4)        % y'
    2+u(2)];      % y'' = 2+x'
[t,u] = ode45(f,tspan,u0);
plot(t,u)

[Devansh Gupta]

Thanx.

If we have an equation of the form

with same boundary conditions then can we use the ode45 function to solve the equation?

[darova]

Is known?

[Devansh Gupta]

It could be any function of x but in this case let it be a number.

[darova]

Yes, the answer is yes. You need to express and use ode45

[Devansh Gupta]

I've expressed it like this

cons1 = ((A(1,1)*D(1,1))-(B(1,1)^2))/(A(1,1));
cons2 = (B(1,1)/A(1,1));
syms w(x)
Dq = diff((q*x),x);
Dw = diff(w,x);
D2w = diff(w,x,2);
D3w = diff(w,x,3);
D4w = diff(w,x,4);
ode = cons1 * D4w == -cons2 * Dq;

but I'm unable to apply the boundary conditions to this expression.

Answer 1

[David Goodmanson]

Hello Devansh

Since you have boundary conditions at both x=0 and x=L you can use bvp4c to solve this. But I believe you have too many boundary conditions. Dropping the subscrips on A,B,D and denoting px by p and pz by q, then

Au''  - Bw'''   + p(x) = 0          (1)
Bu''' - Dw''''  + q(x) = 0          (2)

and substituting as before gives

w'''' = (Aq(x) - Bp'(x))/(AD-B^2)     (3)
u'' =   (Bw'''-p(x))/A                (4)  

( u''' is not involved any more since the substitution process assures that if (3) and (4) work, so does (2) ).

There is a fourth order and a second order equation, so there will be six bc's and not eight. If you keep

w(0) = w'(0) = w(L) = w'(L) = 0,

then there are two bc's left for u. The code below assumes

u(0) = u(L) = 0

but see the note at the end. I assumed p, p' and q to be known algebraic functions.

L = 1.3;
% p,q,A,B,D are in function below
xvec = linspace(0,L,100);
solinit = bvpinit(xvec, @guess);
sol = bvp4c(@fun, @bcs, solinit);
x = sol.x';
% wu ~~ [w w' w'' w''' u u']
wu = sol.y';
w = wu(:,1:4);
u = wu(:,5:6);
figure(1); grid on
plot(x,wu)
legend('w','w''','w''''','w''''''','u','u''')

function dbydx = fun(x,wu)
% wu ~~ [w w' w'' w''' u u']
A = 1;
B = 2;
D = 1;
q = cos(7*x);
p = x^(3/2);
pprime = (3/2)*x^(1/2);
wiv = (A*q-B*pprime)/(A*D-B^2);
uii = (B*wu(4) -p)/A;
dbydx = [wu(2:4); wiv; wu(6); uii];
end
function bcvec = bcs(ya,yb)
bcvec = [ya([1 2 5]); yb([1 2 5])];
end
function g = guess(x)
g = [0 0 0 0 0 0]';
end

For bc's on u, you can specify values of u at both ends, or a value of u at one end and a value of u' at the other, but you can't specify values of u' at both ends. That's because integrating both sides of (4) from 0 to L gives

u'(L) - u'(0)  =  integral(right hand side)

and that condition can't be met with arbitrary choices of u' at each end.