Solving boundary value problem for piecewise defined differential equation

Question

Evgeniy 2014-1-14

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/112314-solving-boundary-value-problem-for-piecewise-defined-differential-equation

回答： Evgeniy 2014-2-12

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.

%

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Evgeniy 2014-2-12

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/112314-solving-boundary-value-problem-for-piecewise-defined-differential-equation#answer_124043

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Answer 2

Mischa Kim 2014-1-14

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/112314-solving-boundary-value-problem-for-piecewise-defined-differential-equation#answer_120823

编辑：Mischa Kim 2014-1-14

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

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.
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 .
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 .

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

Mischa Kim 2014-1-16

编辑：Mischa Kim 2014-1-16

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 2014-1-17

在 MATLAB Online 中打开

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 2014-1-15

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/112314-solving-boundary-value-problem-for-piecewise-defined-differential-equation#answer_120903

编辑：Evgeniy 2014-1-15

在 MATLAB Online 中打开

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,:));

%

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Answer 4

Evgeniy 2014-1-15

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/112314-solving-boundary-value-problem-for-piecewise-defined-differential-equation#answer_120908

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.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Answer 5

Bjorn Gustavsson 2014-1-17

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/112314-solving-boundary-value-problem-for-piecewise-defined-differential-equation#answer_121176

在 MATLAB Online 中打开

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)