MATLAB Answers

# solving simultaneous equations of fourth order PDE in MATLAB

8 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?
##### 10 CommentsShowHide 9 older comments
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.

Sign in to comment.

### Answers (1)

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.
##### 0 CommentsShowHide -1 older comments

Sign in to comment.

R2019a

### Community Treasure Hunt

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

Start Hunting!