Two Nonlinear Coupled PDE's with Neumann BC
2 次查看(过去 30 天)
显示 更早的评论
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 个评论
Jitendra Kumar Singh
2017-10-27
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
2017-10-27
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
2017-10-27
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
2017-10-27
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
2017-10-27
@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
2017-10-27
编辑:Torsten
2017-10-27
@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
2017-10-27
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
2017-10-27
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
2017-10-27
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
2017-10-28
编辑:Jitendra Kumar Singh
2017-10-28
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.
回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 PDE Solvers 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!发生错误
由于页面发生更改,无法完成操作。请重新加载页面以查看其更新后的状态。
您也可以从以下列表中选择网站:
如何获得最佳网站性能
选择中国网站(中文或英文)以获得最佳网站性能。其他 MathWorks 国家/地区网站并未针对您所在位置的访问进行优化。
美洲
- América Latina (Español)
- Canada (English)
- United States (English)
欧洲
- 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)
亚太
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)