# solving simultaneous equations of fourth order PDE in MATLAB

67 views (last 30 days)
Devansh Gupta on 13 Mar 2021
Edited: David Goodmanson on 18 Mar 2021
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?
Devansh Gupta on 15 Mar 2021
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.

David Goodmanson on 17 Mar 2021
Edited: David Goodmanson on 18 Mar 2021
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.

R2019a

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!