You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
Two Nonlinear Coupled PDE's with Neumann BC
1 view (last 30 days)
Show older comments
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.
12 Comments
Jitendra Kumar Singh
on 27 Oct 2017
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;
Jitendra Kumar Singh
on 27 Oct 2017
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
on 27 Oct 2017
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
on 27 Oct 2017
@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
on 27 Oct 2017
Edited: Torsten
on 27 Oct 2017
@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
on 27 Oct 2017
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);
Jitendra Kumar Singh
on 27 Oct 2017
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
on 27 Oct 2017
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
on 28 Oct 2017
Edited: Jitendra Kumar Singh
on 28 Oct 2017
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.
Answers (0)
See Also
Categories
Find more on General PDEs in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)