System of 2 PDEs

2 views (last 30 days)
Olivier L.
Olivier L. on 25 Aug 2015
Commented: Torsten on 27 Aug 2015
Hello,
I need to solve this PDE:
[m.(∂^2 u)/∂t^2 +B.∂u/∂t] . √(1+(∂u/∂x)^2 ) = T.∂/∂x [(∂u/∂x)/√(1+(∂u/∂x)^2 )] + b.M.σ_ext
With m, B, T, b, M and σ_ext constants.
My guess is I need to turn this equation into a system of 2 first-order PDEs:
1) (∂u_1)/∂t=u_2
2) m.(∂u_2)/∂t . √(1+(∂u_1/∂x)^2 ) = T.∂/∂x [(∂u_1/∂x)/√(1+(∂u_1/∂x)^2 )] + b.M.σ_ext-Bu_2.√(1+(∂u_1/∂x)^2 )
Using the format required by pdepe, I rewrite these equations with the following factors c, b and s:
1) C_1=1 ; b_1=0 ; S_1=u_2
2) C_2= m.√(1+(∂u_1/∂x)^2 ) ; b_2= T.(∂u_1/∂x)/√(1+(∂u_1/∂x)^2 ) ; S_2=b.M.σ_ext-Bu_2 √(1+(∂u_1/∂x)^2 )
Can anyone please confirm that this is the correct method?
As for the initial conditions, the 2 functions must be equal to 0 over the space interval x at t=0: u_1(x, t=0) = 0 u_2(x, t=0) = 0
As for the boundary conditions, u_1 and u_2 must be equal to 0 at any t, on the left and right limits, so I wrote p1=q1=0 and p2=q2=0. I do have a doubt about that though.
When I solve the problem, I get the error message "This DAE appears to be of index > 1 "
Can anyone please help me?
Thanks in advance, Olivier
  1 Comment
Torsten
Torsten on 25 Aug 2015
At least the boundary conditions are not implemented correctly.
If you want to set u_1 = u_2 = 0 at both ends,
pl=[ul(1) ; ul(2)]
ql=[0 ; 0];
pr=[ur(1) ; ur(2)];
qr=[0 ; 0];
If your PDE has a name, I'd search whether there are some aspects that should be taken into account for the numerical solution technique.
Best wishes
Torsten.

Sign in to comment.

Answers (1)

Olivier L.
Olivier L. on 26 Aug 2015
Thank you Torsten for your answer regarding the boundary conditions.
I can run my code now, however u_1(x) and u_2(x) remain desperately equal to 0 at any given time and I do not understand why. Is my code incorrect?
Here are my 4 files:
1) [bc2.m]
function [pl,ql,pr,qr] = bc2(xl,ul,xr,ur,t)
pl = [ul(1);ul(2)];
ql = [0;0];
pr = [ur(1);ur(2)];
qr = [0;0];
2) [eqn2.m]
function [c,b,s] = eqn2(x,t,u,DuDx)
mm=1;
Ttheta = 1;
b = 1;
Mp = 1;
sigmaext = 1;
B = 1;
c = [1; mm*(1+DuDx(1).^2).^0.5];
b = [0; Ttheta/(1+DuDx(1).^2).^0.5 * DuDx(1)];
s = [u(2); b*Mp*sigmaext-B*u(2)*(1+DuDx(1).^2).^0.5];
3) [initial2.m]
function value = initial2(x);
value = [0; 0];
4) [eqn2.m]
m = 0;
x = linspace(0,100,100);
t = linspace(0,10,10);
sol = pdepe(m,@eqn2,@initial2,@bc2,x,t);
u1 = sol(:,:,1);
u2 = sol(:,:,2);
surf(x,t,u1);
title('u1(x,t)');
Thanks in advance for your help!
Olivier
  2 Comments
Torsten
Torsten on 26 Aug 2015
I don't see an error in the implementation.
Maybe you could test what happens if you take a value bigger than 1 for b*Mp*sigmaext.
Best wishes
Torsten.
Torsten
Torsten on 27 Aug 2015
I guess the boundary conditions you impose are insufficient to fix a solution.
What you impose at both ends of the interval is
u1(t)=0 and du1(t)/dt = 0.
But if you set u1(t)=0 for all times t, du1(t)/dt = 0 is satisfied automatically. Thus the problem looks underdetermined for me.
I repeat that you should consult references about the numerical solution of this complicated PDE.
Best wishes
Torsten.

Sign in to comment.

Tags

Community Treasure Hunt

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

Start Hunting!