Second order ODE - BVP

Question

sko 2021-2-1

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/732968-second-order-ode-bvp

评论： sko 2021-2-2

采纳的回答： Alan Stevens

在 MATLAB Online 中打开

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,

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

请先登录，再进行评论。

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

Answer 1

Alan Stevens 2021-2-1

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/732968-second-order-ode-bvp#answer_612068

在 MATLAB Online 中打开

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

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

Alan Stevens 2021-2-2

编辑：Alan Stevens 2021-2-2

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 2021-2-2

OMG

:)

Nothing else to add!

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

请先登录，再进行评论。

Answer 2

Alan Stevens 2021-2-1

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/732968-second-order-ode-bvp#answer_612008

编辑：Alan Stevens 2021-2-1

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.

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

sko 2021-2-1

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

请先登录，再进行评论。