0 个投票

Hello. I'm trying to plot a differntial equation using MATLAB. When I run my program, I'm getting different graph.
I've attached the differential equation and my MATLAB coding below.
Can anyone tell me how to use the boundary conditions correctly?
function abcd
ex4init=bvpinit(linspace(0.00001,100,21),[0,0]);
sol=bvp4c(@ex4ode,@ex4bc,ex4init)
plot(sol.x,sol.y(1,:),'blue');
end
function f=ex4ode(r,y)
a=0.001391304348;
b=19.31459475;
f=[y(2)
-(2/r)*y(2)-a*(1-exp(b*y(1)))
];
end
function res=ex4bc(ya,yb)
res=[
ya(1)+1
yb(2)
]
end
This is the graph which I got.
But this is the exact graph.

2 个评论
显示 无 隐藏 无

Robert U
Robert U 2020-3-4
Hi Joy Salomi,
Check your x axis range. I suppose the values should be in nm, i.e. 1e-9 m.
Kind regards,
Robert
Joy Salomi
Joy Salomi 2020-3-4
Hello Sir,
If I change x axis range to 1e-9 m, I didn't get any graph. Plot is empty.
What shall I do?

请先登录，再进行评论。

请先登录，再回答此问题。

 采纳的回答

Robert U
Robert U 2020-3-5
编辑：Robert U 2020-3-6

0 个投票

Hi Joy Salomi,
Thanks for the paper. The differential equation should be written in spherical coordinates (edit: I think there is a mistake to assume R to be constant). Since R is assumed to be constant there is no singularity. I adjusted the code and tested it in Matlab 2019b. Depending on chosen xmesh results might differ. You have to try a bit to get a stable solution. I attached an example for the D = 200 nm graph.
abcd.m
function abcd(ah,D,phi_s)
ex4init=bvpinit(linspace(0,D/2,3e3),[phi_s 0]);
sol=bvp4c(@ex4ode,@ex4bc,ex4init);
plot(ah,sol.x*1e9,sol.y(1,:));
ah.XLabel.String = 'r [nm]';
ah.YLabel.String = 'phi [V]';
function dphidr=ex4ode(r,phi)
a = -1.602176634e-19*1e17*(1e-2)^(-3)/(8.8541878125e-12*11.5);
b = 1.602176634e-19/(physconst('Boltzmann')*600);
dphidr = [ phi(2)
a * ( 1-exp( b*phi(1) ) ) ];
end
function res=ex4bc(ya,yb)
res=[ ya(1)+ phi_s
yb(2) ];
end
end
Example call:
fh = figure;
ah = axes(fh);
hold(ah,'on')
arrayfun(@(dIn) abcd(ah,200e-9,dIn),[0.1,0.3,0.7,1])
Result:
Kind regards,
Robert

4 个评论
显示 2更早的评论 隐藏 2更早的评论

Joy Salomi
Joy Salomi 2020-3-6
Hello Sir,
Thank you for your code.
What is 'ah' in the above code?
R is not constant Sir.
That equation (23) can be written as phi''(r)+(2/r)phi'(r)=-a(1-e^(b*phi(r)));
This is the equation sir.
Robert U
Robert U 2020-3-6
编辑：Robert U 2020-3-6
Hi Joy Salomi,
'ah' is the axis handle to plot the graph into. That helps to draw multiple lines into one plot.
If you run my code you can easily see, that the figure of the given paper is reproduced. Changing D to 50e-9 leads to the 50 nm grain size figure.
The variable 'R' is not defined within the paper (sic). Since 'r' and 'R' will not be the same value, I assumed R = D/2. That does not correspond to the Poisson-Boltzmann-equation formulated here. But I am quite sure it is what has been solved.
Kind regards,
Robert
Robert U
Robert U 2020-3-6
Hi Joy Salomi,
I added the singular term as you proposed according to my first answer:
abcd_wSingularity:
function abcd_wSingularity(ah,D,phi_s)
ex4init=bvpinit(linspace(0,D/2,3e3),[phi_s 0]);
S = [0 0; 0 -2];
options = bvpset('SingularTerm',S);
sol=bvp4c(@ex4ode,@ex4bc,ex4init,options);
plot(ah,sol.x*1e9,sol.y(1,:));
ah.XLabel.String = 'r [nm]';
ah.YLabel.String = 'phi [V]';
function dphidr=ex4ode(r,phi)
a = -1.602176634e-19*1e17*(1e-2)^(-3)/(8.8541878125e-12*11.5);
b = 1.602176634e-19/(physconst('Boltzmann')*600);
dphidr = [ phi(2)
a * ( 1-exp( b*phi(1) ) ) ];
end
function res=ex4bc(ya,yb)
res=[ ya(1)+ phi_s
yb(2) ];
end
end
Example call:
fh = figure;
ah = axes(fh);
hold(ah,'on')
arrayfun(@(dIn) abcd_wSingularity(ah,200e-9,dIn),[0.1,0.3,0.7,1])
Result:
Kind regards,
Robert
Joy Salomi
Joy Salomi 2020-3-6
编辑：Joy Salomi 2020-3-6
Hello Sir,
I understood now. Highly satisfied with your answer.
The code you sent is working.
Thank you so much for your help Sir. :)

请先登录，再进行评论。

更多回答（2 个）

0 个投票

This is a bit hard to solve without knowing more of the problem. What I would do is to break down the problem and analyse each step. For example, look at the values you have in the solution "sol". I plot separately:
>> plot(sol.x);
>> plot(sol.y);
And also
>> plot(sol.y');
Is anything here remotely connected to what you are expecting? My guess is that no, sol.y has values between 0 and 12e4 whilst the graph you want is between -1 and 0. So the problem is not really about Matlab but about the equations that you are using.
Hope this helps.

1 个评论
显示 -1更早的评论 隐藏 -1更早的评论

Joy Salomi
Joy Salomi 2020-3-4
编辑：Joy Salomi 2020-3-4
Thank you for your reply Sir. But my doubt is somewhat different. I'm trying to get the exact curve which is already available.
This is the curve I want.

请先登录，再进行评论。

Robert U
Robert U 2020-3-4

0 个投票

Hi Joy Salomi,
the implementation of a boundary value problem with singular term has to be done a bit different than you did.
--> doc: Solve BVP with Singular Term
function abcd
ex4init=bvpinit(linspace(0,100,5),[-1 0]);
S = [0 0; 0 -2];
options = bvpset('SingularTerm',S);
sol=bvp4c(@ex4ode,@ex4bc,ex4init,options);
plot(sol.x,sol.y(1,:),'blue');
end
function f=ex4ode(r,y)
a=0.001391304348;
b=19.31459475;
f = [ y(2)
-a*(1-exp(b*y(1))) ];
end
function res=ex4bc(ya,yb)
res=[ ya(1)+1
yb(2) ];
end
Using the values "a" and "b" you supplied still leads to a singular Jacobian matrix which provoces an error. Changing both these values to one leads to a solution that qualitatively corresponds to the curves drawn.
There are two more points that hinder me to test my solution:
  1. I don't know how "a" and "b" translate to "phi_s".
  2. I still assume that r is given in SI units. If r is in unit [m] the interval needs to be adjusted to [nm] as shown in your graph.
Kind regards,
Robert

4 个评论
显示 2更早的评论 隐藏 2更早的评论

Joy Salomi
Joy Salomi 2020-3-5
编辑：Joy Salomi 2020-3-5
Hello Sir,
Thank you for your reply.
One of the boundary conditions is y(0)= - phi_s
I have taken phi_s as 1 and mentioned that as y(0)=-1 in my problem.
Joy Salomi
Joy Salomi 2020-3-5
编辑：Joy Salomi 2020-3-5
Hello Sir,
Can I use the following code? Still I didn't get the exact curve. Is my way of using boundary conditions in the below coding correct?
Boundary conditions are: y(0)=-1 , y'(100)=0
function examplesperical
m =2;
x =linspace(0,25);
t=linspace(0,100);
sol = pdepe(m,@pdex4pde,@pdex4ic,@pdex4bc,x,t);
u1 = sol(:,:,1);
figure
plot(x,u1(end,:))
%------------------------------------------------------------------
function [c,f,s] = pdex4pde(x,t,u,DuDx)
c = 1;
f = DuDx;
a=0.001391304348;
b= 19.31459475;
F =-a*(1-exp(b*u(1)));
s = F;
% --------------------------------------------------------------
function u0 = pdex4ic(x)
%create a initial conditions
u0 = -1;
% --------------------------------------------------------------
function[pl,ql,pr,qr]=pdex4bc(xl,ul,xr,ur,t) %create a boundary conditions
pl = ul(1)+1;
ql = 0;
pr = ur(1)-0;
qr = 1;
Robert U
Robert U 2020-3-5
Hi Joy Salomi,
maybe you can state what physical process is described by that differential equation and what reference you are using (where is your photo from).
I still think that your values for a and b do not correspond to the case that is depicted in the book/article.
Kind regards,
Robert
Joy Salomi
Joy Salomi 2020-3-5
编辑：Joy Salomi 2020-3-5
Hello Sir,
I have attached the PDF for your reference.
a and b are combination of some parametrs.
By subtituting the values of the parameters, I have found the values of a and b.

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Programming 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by