Solving two second order BVPs and the boundary conditions are related

Question

Mohammed Elhefnawy 2019-12-13

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/496537-solving-two-second-order-bvps-and-the-boundary-conditions-are-related

评论： Mohammed Elhefnawy 2019-12-15

采纳的回答： darova

在 MATLAB Online 中打开

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 个评论
显示 2更早的评论隐藏 2更早的评论

darova 2019-12-14

Solve for region 2 (-b .. 0)

Get solution and solve for region 1 (0 .. a)

Mohammed Elhefnawy 2019-12-14

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

请先登录，再进行评论。

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

Answer 1

darova 2019-12-14

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/496537-solving-two-second-order-bvps-and-the-boundary-conditions-are-related#answer_406362

在 MATLAB Online 中打开

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