Solving boundary value problem for piecewise defined differential equation

Greetings. I faced the problem with solving the system of differential equations showed below. This system contains two equations (one for each interval). Due to the fact that these are 2-nd order differential equations and I don't know the values of derivatives at the ends (y'(-20) and y'(10)) this is boundary problem and the system could be solved by bvp4c function.

The problem is how to define the condition of the first derivatives near the 0 (bottom of the picture). I tried to do this:

function [ dydx ] = tode( x , y )
dydx = zeros(2,1); % create zero array
dydx(1) = y(2);
if (x > 0)
    dydx(2) = b/a*((1+y(1))^(3/2)-1);
elseif (x <= 0)
    dydx(2) = 0;
end;
end

But this expression doesn't include the condition of first-order derivatives. How can I solve the problem? Thank you.

%

Answer 1

[Evgeniy]

The both bvp4c and bpv5c solvers can find the solution for piecewise defined equations. In this case the solution must be smooth at each region. At the borders you can specify additional conditions. Here is the link which describes how to do that: http://www.mathworks.com.au/help/matlab/math/boundary-value-problems.html

Answer 2

[Mischa Kim]

Look at this problem as an initial value problem and turn it into an optimization problem:

1. Set up an algorithm that integrates the piece-wise function from the left, that is in [-20,0], and, separately, from the right, in [10,0]. The only unknowns are the two first derivatives of y at the boundaries (let's call them yp = [yp1, yp2] ), correct? So you have to make a guess.
2. Once you have performed the integration, the two boundary conditions at x=0 (most likely) will not be satisfied: dy(0)≠ 0, dyp(0)≠ 0 .
3. You should then be able to use a root-finding algorithm (e.g. Newton-Raphson) to find the root of the cost function J of the optimization, which in this case is J = norm([dy(0);dyp(0)]). Note, when both boundary conditions at x=0 are satisfied, J=0 .

[Evgeniy]

I can say with 99% of confidence that y'(10)=0. In this case I should change the first derivative at the boundary x=-20 until it satisfy the two described conditions. But in this case the solution of the positive part will not be affected by changing of initial values of negative part (if for example y(-20)=3) but it should. How should I avoid this? Thank you.

PS And also, what did you mean by 'dyp'?

[Mischa Kim]

Well, in this case the problem becomes much easier. Solve the initial value problem in x = [10,0] (using, e.g., ode45). This will give you y(0+) and y'(0+), which you can use for the initial conditions for the second integration in x = [0,-20]. Problem solved.

In fact, since the 2nd differential equation has such an easy form, you can do this part by hand.

Note, however, that y(-20) = 5, which needs to be confirmed by the second integration. If the result differs your previous assumption (99%) was not correct. Are all the constants known?

[Evgeniy]

The value of function at point 10 must be equal to 0 (y(10)=0) because it is potential of the metal. Also the first derivative at this point should be equal to zero or near that value. I tried to solve these equations using the method described above. The code is this:

global coef_A coef_B;

coef_A = 1.5477e+10;
coef_B = 9.7420e+03;

U = 1;              %value at y(-20e-7)
dy1 = 0;           %value at y'(-20e7) (a guess)

err = 1e-1;       %allowed error for value at y(10e-7)
end_val = 2;     %value at y(10e-7)

xspan_neg = [-20e-7 : 1e-8 : -1e-8];      %defining intervals
xspan_pos = [1e-8 : 1e-8 : 10e-7];

while end_val > err         %while and value is bigger than the error

      y0_neg = [U dy1];      %initial values at -20e-7
      [X1,Y1] = ode45('neg_diff_eq',xspan_neg,y0_neg);  %finding the y(x) and y'(x) at interval [-20e-7; 0]

      U1 = Y1(end,1);          %value at y(+0), first continuity condition
      dU1 = 0.5391*Y1(end,2);  %value at y'(+0), second continuity condition 

      y0_pos = [U1 dU1];       %initial values at y(+0) and y'(+0)
      [X2,Y2] = ode45('pos_diff_eq',xspan_pos,y0_pos);   %finding the y(x) and y'(x) at interval [0; 10e-7]

      end_val = abs(Y2(end,1));   %found value at y(10e-7)
      dy1 = dy1 - 1e4;            %decreasing the value of first derivative at point -20e-7

end

X = [X1; X2]; Y = [Y1; Y2];      %combining results

subplot(2,1,1);                  %drawing
plot(X,Y(:,1));
subplot(2,1,2);
plot(X,Y(:,2));

And functions are look like:

function [ dydx ] = neg_diff_eq( x , y )

dydx = zeros(2,1); % create zero array
dydx(1) = y(2);
dydx(2) = 0;

end

function [ dydx ] = pos_diff_eq( x , y )

global coef_A coef_B;

dydx = zeros(2,1); % create zero array
dydx(1) = y(2);

dydx(2) = coef_A*(1+coef_B*y(1))^(3/2)-1;

end

But when I run the computation MATLAB just freezes and shows 'Busy' in the bottom of the main window. What is wrong with this code?

Answer 3

[Evgeniy]

Thanks for the answer. This problem is about the finding the potential inside the capacitor. The boundary conditions y(-20)=5 and y(10)=0 are potentials of capacitor plates. The capacitor is filled with dielectric (x = [-20;0]) and semiconductor (x = [0;10]). The potential inside the dielectric should linearly fall from 5V to y(-0). Inside the semiconductor due to accumulation of charges near the boundary y(+0) the potential should rise while being almost zero elsewhere (see picture below). There are conditions of continuity of potential (y(-0)=y(+0)) and electric field displacement (a*y'(+0)=c*y'(-0)).

I have tried to solve these two equations separately but I didn't get the success. The thing is that I don't know the first derivative (electric field) of both dielectric and semiconductor region. Regarding the dielectric the first derivative should be equal to y'(x)=(a/c)*y'(+0)=const.

I also tried to use bvp4c function and got the result very similar to the graph of capacitor but the derivative near 0 has very strange behavior (see picture). I think that is because I didn't specify the condition of electric field displacement continuity. The Matlab code is:

function [ dydx ] = tode( x , y )
dydx = zeros(2,1); % create zero array
dydx(1) = y(2);
if (x > 0)
    dydx(2) = b/a*((1+y(1))^(3/2)-1);
elseif (x <= 0)
    dydx(2) = 0;
end;
end

function [ res ] = tobc( ya , yb )
U1 = 5;
U2 = 0;

res = [ ya(1)-U1
        yb(1)-U2];
end

solinit = bvpinit(linspace(-20,10,300),[5 0]);

sol = bvp4c(@tode, @tobc, solinit);

x = linspace(-20,10,300);
y = deval(sol,x);
plot(x,y(1,:));

%

Answer 4

[Evgeniy]

Maybe because the bvp4c may be applicable only if the function and derivative are smooth but due to condition a*y'(+0)=c*y'(-0) it seems that first derivative is piecewise smooth so it is impossible to use bvp4c to solve this equation.

Answer 5

[Bjorn Gustavsson]

Second differential equation should give you a simple first order polynomial:

y = 5 + C * (20 + x)

That should give you simple expressions of both y(-0) and dy/dx(-0). That should be possible to solve with the suggested shooting method for only the possitive ODE - but taking into account both y(+0) and dy/dx(+0)