Two Nonlinear Coupled PDE's with Neumann BC

Does MATLAB provide a solver for Nonlinear Coupeled PDE's that I have attached in the image file. If yes, then which solver can be used. As I am new to MATLAb, can anyone provide a link to same sort of example.

The equations are in the image file.

The research paper if the model is:- http://aip.scitation.org/doi/abs/10.1063/1.1316047

[Torsten]

Check out "pdepe".

What is the domain for z ?

Best wishes

Torsten.

[Jitendra Kumar Singh]

Thank you for your reply.

z ranges from 0 to infinity. At infinity value of N and T are zero. However I want to model only for finite z.

I tried * pdepe*, but it gave me some unusual error. See my code for first PDE and the associated error.

*??? Error using ==> daeic12 at 77 This DAE appears to be of index greater than 1.

Error in ==> ode15s at 395 [y,yp,f0,dfdy,nFE,nPD,Jfac] = daeic12(odeFcn,odeArgs,t,ICtype,Mt,y,yp0,f0,...

Error in ==> pdepe at 320 [t,y] = ode15s(@pdeodes,t,y0,opts);

Error in ==> trial at 7 sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);*

 function pdex1

 m = 0;
 x = linspace(0,1,20);
 t = linspace(0,2,5);

 sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);
 % Extract the first solution component as u.
 N = sol(:,:,1);

 % A surface plot is often a good way to study a solution.
surf(x,t,N) 
title('Numerical solution computed with 20 mesh points.')
xlabel('Distance x')
ylabel('Time t')

 % --------------------------------------------------------------
function [c,f,s] = pdex1pde(x,t,N,dNdx)
c = 1;
D0 = 18;
A = 7*10^18;
B = 0.54*10^18;
D = D0*A/(A + D0*N*abs(log(1+(B/N)^(2/3))));
f = D.*dNdx;
s = 0;
% --------------------------------------------------------------
function N0 = pdex1ic(x)
R0 = 0.33;
fp = 1.25;
hv = 1.55;
delta = 9.8;
N0 = ((1-R0)*fp/(hv*delta))*exp(-x/delta);

 % --------------------------------------------------------------
function [pl,ql,pr,qr] = pdex1bc(xl,Nl,xr,Nr,t)
R0 = 0.33;
fp = 1.25;
hv = 1.55;
delta = 9.8;

 S = 3*10^5;

 D0 = 18;
A = 7*10^18;
B = 0.54*10^18;
D = D0*A/(A + D0*Nl*abs(log(1+(B/Nl)^(2/3))));

 pl = -((S*fp*(1-R0))/(hv*delta))*((1/D0) + (Nl/A)*abs(log(1+(B/Nl)^(2/3))));
ql = 0;
pr = Nr;
qr = 1;

[Torsten]

Please write out the boundary conditions you want to specify in a mathematical notation so that it is possible to compare with what you actually set.

Best wishes

Torsten.

[Jitendra Kumar Singh]

As at right boundary I do not know the value of N, so I have set it as Nr because MATLAB may approximate it itself. However on the left boundary according to MATLAB format(value+derivative|@boundary = 0; pl or pr = value), I am using pl= -S*N(0,t)*[D(N)^(-1)]. Also as I don't know the value of N(0,t), so I placed value of N(0,0) where it should be N(0,t). This might be wrong but the solver should have given some result. Although the result can be weird.

[Torsten]

As written in your function file, you set

@x=0:  -((S*fp*(1-R0))/(hv*delta))*((1/D0) + (N/A)*abs(log(1+(B/N)^(2/3)))) = 0

@x=20: N + (D0*A/(A + D0*N*abs(log(1+(B/N)^(2/3)))))*dN/dx = 0

Is this really what you want to set ?

Best wishes

Torsten.

[Jitendra Kumar Singh]

@x=20: I don't know the value, so I want MATLAB to pick an approximate value by itself.

@x=0: I want the value to be +((S*fp*(1-R0))/(hv*delta))*((1/D0) + (N/A)*abs(log(1+(B/N)^(2/3)))) There is no "equal to zero" relation there. According to matlab format, The equation @x=0 would read,

+[((S*fp*(1-R0))/(hv*delta))*((1/D0) + (N/A)*abs(log(1+(B/N)^(2/3))))] + {[dN/dx]@x=0} = 0

[Torsten]

@x=20: Matlab can't "pick" a value by itself. If you do not know the value for N, then setting dN/dx = 0 often is an option:

pr=0.0, qr=1.0

@z=0: You set

pl = -((S*fp*(1-R0))/(hv*delta))*((1/D0) + (Nl/A)*abs(log(1+(B/Nl)^(2/3))))
ql = 0

The boundary condition is

p+q*f = 0

Since

f=(D0*A/(A + D0*N*abs(log(1+(B/N)^(2/3)))))*dN/dx,

you set

-((S*fp*(1-R0))/(hv*delta))*((1/D0) + (N/A)*abs(log(1+(B/N)^(2/3)))) + 0*(D0*A/(A + D0*N*abs(log(1+(B/N)^(2/3)))))*dN/dx)=0

thus

-((S*fp*(1-R0))/(hv*delta))*((1/D0) + (N/A)*abs(log(1+(B/N)^(2/3)))) = 0.

In the example,

f=du/dx
pr=pi*exp(-t)
qr=1.0,

thus

p+q*f = pi*exp(-t)+1*du/dx = 0.

Best wishes

Torsten.

[Jitendra Kumar Singh]

Thank you so much for your time and help. You indeed made many things clear to me. Just thank you wouldn't be enough.

However, the errors are still there. Does it have to do with the thing that variable N comes in the boundary conditions and pdepe isn't made for such BC's.

The errors are:-

 ??? Error using ==> daeic12 at 77
This DAE appears to be of index greater than 1.
 Error in ==> ode15s at 395
    [y,yp,f0,dfdy,nFE,nPD,Jfac] = daeic12(odeFcn,odeArgs,t,ICtype,Mt,y,yp0,f0,...
 Error in ==> pdepe at 320
  [t,y] = ode15s(@pdeodes,t,y0,opts);
 Error in ==> trial at 7
sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);

[Torsten]

Please show the updated code.

Best wishes

Torsten.

[Jitendra Kumar Singh]

 function pdex1

 m = 0;
x = linspace(0,1,20);
t = linspace(0,2,5);

 sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);
% Extract the first solution component as u.
N = sol(:,:,1);

 % A surface plot is often a good way to study a solution.
surf(x,t,N) 
title('Numerical solution computed with 20 mesh points.')
xlabel('Distance x')
ylabel('Time t')

 % --------------------------------------------------------------
function [c,f,s] = pdex1pde(x,t,N,dNdx)
c = 1;
D0 = 18;
A = 7*10^18;
B = 0.54*10^18;
D = D0*A/(A + D0*N*abs(log(1+(B/N)^(2/3))));
f = D.*dNdx;
s = 0;
 % --------------------------------------------------------------
function N0 = pdex1ic(x)
R0 = 0.33;
fp = 1.25;
hv = 1.55;
delta = 9.8;
N0 = ((1-R0)*fp/(hv*delta))*exp(-x/delta);

 % --------------------------------------------------------------
function [pl,ql,pr,qr] = pdex1bc(xl,Nl,xr,Nr,t)
R0 = 0.33;
fp = 1.25;
hv = 1.55;
delta = 9.8;

 S = 3*10^5;

 D0 = 18;
A = 7*10^18;
B = 0.54*10^18;

 pl = -((S*fp*(1-R0))/(hv*delta))*((1/D0) + (Nl/A)*abs(log(1+(B/Nl)^(2/3))));
ql = 0;
pr = 0;
qr = 1;

[Torsten]

What approxiamte value for N @z=0 do you get from

-((S*fp*(1-R0))/(hv*delta))*((1/D0) + (Nl/A)*abs(log(1+(B/Nl)^(2/3))))=0

?

Is it reasonable ?

Better prescribe this value directly (pl=Nl-value, ql=0), not in an implicit formulation.

Also start with a D which does not depend on N.

Best wishes

Torsten.

[Jitendra Kumar Singh]

Hello, I realised there were many things wrong in my code. The worst were the badly scaled values. Finally, with your help I managed to solve it. I couldn't ever thank you enough for your help.

But still, thanks a ton.