Solving two second order BVPs and the boundary conditions are related
2 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Mohammed Elhefnawy
am 13 Dez. 2019
Kommentiert: Mohammed Elhefnawy
am 15 Dez. 2019
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 Kommentare
Akzeptierte Antwort
darova
am 14 Dez. 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 Kommentare
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Eigenvalue Problems finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!