solving simultaneous equations of fourth order PDE in MATLAB

7 views (last 30 days)
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 Comments
darova
darova on 15 Mar 2021
Yes, the answer is yes. You need to express and use ode45
Devansh Gupta
Devansh Gupta on 15 Mar 2021
Edited: 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
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.

Products


Release

R2019a

Community Treasure Hunt

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

Start Hunting!