MATLAB Answers

Solving two second order BVPs and the boundary conditions are related

2 views (last 30 days)
I am trying to sove these two BVP equations (solving for time independent schrodinger equation in 1D)
the bounadry conditions are related and
the BCs are
and the analytical solution is (before applying BCs)
I am getting some errors which are shown in the code below
I would appreciate any help please
function dudx = osc(x,u)
h=6.62606896E-34; %% Planck constant [J.s]
hbar=h/(2*pi);
q=1.602176487E-19; %% electron charge [C]
%m=input('enter the mass ');
m=9.10938188E-31; %% electron mass [kg]
%E=input('enter the energy ');
E=10000;
k=(2*m*E)/(hbar^(2));
V=-(1.6*(10^(-19)))/(4*3.14*8.854*(10^(-12))*abs(x));
A=sqrt((2*m*E)/(hbar^(2)));
B=sqrt(((2*m)/(hbar^(2)))*(E-V));
dudx = [u(2)
-2*1i*k*(u(2))+((k^(2))-(A^(2)))*u(1)
u(4) %% error :index exceeds array bounds
2*1i*k(u(4))+((k^(2))-(B^(2)))*u(3)]; %% error: Array indices must be positive integers or logical values.
end
%-----------------------------------------------------------------------%
function res = bcfcn(ya,yb,yc) % boundary conditions
res = [ya(1)-ya(3)
ya(2)-ya(4) %%error : the function bcfcn should return a column of two vectors
yb(1)-yc(1) %% and I am m=not sure if it is right to write the BCs like this
yb(2)-yc(2)];
end
%-----------------------------------------------------------------------%
function g = guess(x) % initial guess for y and y'
g = [exp(x)+exp(-x)
exp(x)+exp(-x)];
end
%-----------------------------------------------------------------------%
xmesh = linspace(-10,10,21);
solinit = bvpinit(xmesh, @guess);
sol = bvp4c(@osc, @bcfcn, solinit);
plot(x,u(:,1));
%-----------------------------------------------------------------------%

  4 Comments

Show 1 older comment
darova
darova on 14 Dec 2019
Solve for region 2 (-b .. 0)
Get solution and solve for region 1 (0 .. a)
Mohammed Elhefnawy
Mohammed Elhefnawy on 14 Dec 2019
Thanks alot
but there is something that I still don'y get it
I will solve for the second region but to solve it I have to put boundary conditions
so if I assume for example that (u2(0))=some constant (A) and the same for the other conditions and solve it
how to take the solution and pass it to the other region and how to deal with the BCs their too

Sign in to comment.

Accepted Answer

darova
darova on 14 Dec 2019
How to set up initial guess? Any ideas?
function main
%-----------------------------------------------------------------------%
b = 10;
a = 10;
h=6.62606896E-34; %% Planck constant [J.s]
hbar=h/(2*pi);
q=1.602176487E-19; %% electron charge [C]
%m=input('enter the mass ');
m=9.10938188E-31; %% electron mass [kg]
%E=input('enter the energy ');
E=10000;
k=(2*m*E)/(hbar^(2));
V = @(x) -(1.6*(10^(-19)))/(4*3.14*8.854*(10^(-12))*abs(x));
A = sqrt((2*m*E)/(hbar^(2)));
B = @(x) sqrt(((2*m)/(hbar^(2)))*(E-V(x)));
myode = @(x,u) [u(2)
-2*1i*k*(u(2))+((k^(2))-(A^(2)))*u(1)
u(4)
-2*1i*k*(u(4))+((k^(2))-(B(x)^(2)))*u(3)];
x0 = fsolve(@F,[1 1 1 1]*exp(-b)); % what initial guess should be set?
[x1,u11] = ode45(myode, [-b a], x0);
plot(x1,u11)
function res = F(x0)
[x,u] = ode45(myode, [-b a], x0);
u1 = u(:,1);
du1 = u(:,2);
u2 = u(:,3);
du2 = u(:,4);
res = [interp1(x,u1,0) - interp1(x,u2,0) % u1(0) == u2(0)
u1(end) - u2(1) % u1(a) == u2(-b)
interp1(x,du1,0) - interp1(x,du2,0) % du1(0)== du2(0)
du1(end) - du2(1)]; % du1(a)== du2(-b)
end
end

  7 Comments

Show 4 older comments
Mohammed Elhefnawy
Mohammed Elhefnawy on 15 Dec 2019
A is (alpha) and B is(Beta)
Annotation 2019-12-15 131345.jpg Annotation 2019-12-15 131350.jpg
Annotation 2019-12-15 131356.jpg
all of this is to try solving time independent schrodinger equation
these equations are form the book (Semiconductor physics and devices by donald A.Neeman)
chapter 3 section 1

Sign in to comment.

More Answers (0)

Sign in to answer this question.

Tags

Products


Release

R2018a

Translated by