Second order ODE - BVP

Hi,

I am trying to obtain the solution of the following second order ODE but I struggle.

and on the interval [0, 0.5]. It is my understanding that 1) transoformation to the system of first order ODE and a Matlab bvp4c solver should be used. I wrote therefore the code:

function bvp5
xlow=0; xhigh=0.5;
solinit = bvpinit(linspace(xlow,xhigh,10),[0 1]); 
sol = bvp4c(@bvp5ode,@bvp5bc,solinit);
xint = linspace(xlow,xhigh);
Sxint = deval(sol,xint);
plot(xint,Sxint(1,:))
% -----------------------------------------------
function dydx = bvp5ode(x,y)
dydx = [ y(2) 0.64*y(1)-(2/x)*y(2) ];
% -----------------------------------------------
function res = bvp5bc(ya,yb)
res = [ ya(1)-0.2 yb(2) ];

I obtain the folelowing error: Unable to solve the collocation equations -- a

singular Jacobian encountered.

1) How to form a guess ? I dont have any idea ...

2) What is the problem with my equation or my code?

Best wishes,

Answer 1

[Alan Stevens]

I used an alternative method that makes use of Matlab's fzero function. You can also see one way of avoiding the problem at x = 0:

xlow=0; xhigh=0.5;
xspan = [xlow xhigh];
dydx0 = 0;

y0 = 0;  % Initial guess for y(x=0)
% Use fzero to get value of y0 that makes y(x=0.5) = 0.1;
y0 = fzero(@(y0) y0val(y0, xspan), y0);
disp(['y(x=0) = ' num2str(y0)])

% Now use ode45 to integrate from x=0 to x=0.5 with good starting value
% for y(x=0)
Y0 = [y0 0]; 
[x, Y] = ode45(@odefn, xspan, Y0);

plot(x,Y(:,1)),grid
xlabel('x'),ylabel('y')

function Z = y0val(y0, xspan)

         Y0 = [y0, 0];
         [~, Y] = ode45(@odefn, xspan, Y0);
         yend = Y(end,1);
         Z = yend - 0.1;
         
end


function dYdx = odefn(x, Y)

         y = Y(1); dydx = Y(2);
         d2ydx2 = 0.64*y;
         if dydx>0
             d2ydx2 = d2ydx2 - (2/x)*dydx;
         end
         dYdx = [dydx;
                 d2ydx2];
         
 end

[sko]

Hey thank you very much. The issue with that approach is that I need to have a good guess which is somewhat an arbitrary choice. I am satisfied with the solution you proposed initially (for a nonzero but very low x). However, I have some inconsistency I d like to discuss. The analytical solution of my ODE exists and when I compare it with the results obtained with bvp4c there is a great disagreement. Can you please help me understand where am I wrong with defining bvp4c problem?

%analytical solution of ODE 
k=6.4;
D=0.1;
% k/D=0.64; 
R=0.5;
y0=0.1;
phi=R*(k/D)^(1/2);
r=linspace(0.01,R,30);
y=y0*(R./r).*(sinh(phi*r/R))./(sinh(phi));
plot(r,y);

As you can see y is being reduced fastly towards the x=0, considering the given ratio of k/D=0.64. Thus, there is some issue with my bvp problem definition ....

NB I d like to point out that shooting method you proposed yields the same results as analytical solution for a proper guess of y. However this method is impractical for my problem.

Thanks

[Alan Stevens]

I think the problem with using bpv4c is that it expects values of y(x=0) and y(x=0.5), but you have y'(x=0), not y(x=0). I'm not very familiar with bvp4c (I rarely use it), so I can't really help you further with it.

It's true that the shooting method requires a reasonable initial guess. Presumably your actual problem is much more complicated if the method is no use to you.

(NB In your analytical solution you have k = 6.4 and D = 0.1. This doesn't result in k/D = 0.64 !!!)

[sko]

OMG

:)

Nothing else to add!

I obtain the good agreement with the analytical solution. Many Many thanks!

Answer 2

[Alan Stevens]

You start at xlow = 0, but your function bvp5ode has an x in the denominator. bvp4c doesn't like it when this is zero! If you set xlow = 0.001, say, then it runs.

[sko]

Hi, Thank you very much for your answer. It shows no error indeed. However, as you may have noticed, my equation is some sort of stationary diffusion-reaction equation, therefore, I should be able to have x=0 (surface of particle). Another important question is the guess function. Could you comment on this subject?

Thank you in advance