Boundary value problem-bvp4c

Hello

I am working on a paper by Di Federico et al 2012, they solved this by Mathematica@. I need to reproduce it in MATLAB and extend it for radial geometry.

Here is the problem :

I am trying to solve an ODE ( has been derived from an PDE by scaling and similarity solution process) which I have two values of it at the end of the interval [0,1] .

I know that as , the function tends to zero .

Also I know that as the derivative equals to .

as "bvp4c" needs the values in two different point, how should I write the residual function "bcfun" for this problem as I do not have two points ?

instead I have information about one end point and the slope near that point

This is the system of first order ODE :

also there is singularity near this point which they handle it by studying the problem in the neighborhood

the solution should have an answer like this :

Answer 1

[Torsten]

Use ode45 and integrate backwards:

https://uk.mathworks.com/matlabcentral/answers/1881197-numerical-integration-backwards-for-ode45

[sepideh]

Thank you

I saw the link. but there is a problem

in the example to sent to me the end points are known (y0=1) (y0=2.7)

in my case I only know that in y(1)=0 and I dont have any information about the starting point

I just know the slop which is y'= -2/3

[Torsten]

Y1 = 0 and Y2 = -2/3 are your initial conditions at x=1 for your system of ODEs. What's the problem ?

And be careful: there is an error in your vector of derivatives: x*Y2 instead of x*Y1 in the numerator of the second component.

[sepideh]

yes you are rigth it is a typo

[Torsten]

Of course - since you divide by Y1 - you cannot set Y1 = 0 at x = 1 exactly. Use some small number instead of 0, e.g. 1e-8.

[sepideh]

Dear Toesten

I wrote 2 codes, non of them work.

Would you please take a look at them and kindly correct them?

[sepideh]

clc

clear

xspan= [0,1]

y0=1

[xf,yf]= ode45(@bvpfun,xspan,y0)

xsapn= [1,0]

y0=1e-10

[xr,yr]= ode45(@bvpfun,xspan,y0)

plot(xf,yf,xr,yr)

function dydx=bvpfun(x,y)

dydx= [y(2)

-(3*y(2)^2+x*y(2)+y(1))/(3*y(1))];

end

[sepideh]

eps= 1e-10;

xmesh= linspace(0,1-eps,10);

solinit= bvpinit(xmesh,[1,0]);

sol=bvp4c(@bvpfun,@bcfun,solinit);

plot(sol.x,sol.y,'o')

function dydx=bvpfun(x,y)

dydx= [y(2),-(3*y(2)^2+x*y(2)+y(1))/(3*y(1))];

end

function res=bcfun(~,yb)

res=[yb(2)+2/3

yb(1)-(2/3*eps)];

end

[Torsten]

xspan = [1 0];
y0 = [1e-8 -2/3];
fun = @(x,y)[y(2);-(3*y(2)^2+x*y(2)+y(1))/(3*y(1))];
[X,Y] = ode45(fun,xspan,y0);
plot(X,Y)

[Torsten]

The blue line is Y1, the orange line is Y2. The matrix Y is made up as Y = [Y1,Y2]. Now you should be able to modify the plot command to only plot Y1.

[sepideh]

thank you

[sepideh]

xspan = [1 0];

y0 = [1e-8 -1/3];

fun = @(x,y)[y(2);-(3*y(2)^2+x*y(2)+y(1))/(3*y(1))];

[X,Y] = ode45(fun,xspan,y0);

plot(X,Y(:,1))

hold on

syms F(Y) %analytical solution

F(Y)= 1/6*(1-Y^2);

fplot(F,[0 1],'o')

axis([0 1 0 1])

hold off

thank you very much

now the analytical and numerical results are perfectly match. you helped me a lot

[Torsten]

Fine. I tried to find the analytical solution with "dsolve", but I didn't succeed. You took it from some publication, or was Mathematica able to find it ?

[sepideh]

I am working on Gravity Currents and the ensuing ODE has a analytical solution for the shape function as : (see Huppert 1982, Di Federico 2012)

this is the case for sudden or instantaneous release (

for the ODE can not solved analytically

here I made a mistake, actually the slope for this situation is -1/3 not -2/3 because now I am trying to solve the ODE for instantaneous release. the code that you gave me was perfectly match with the analytical solution (I correctted the slope) in this situation

when (constant injection) the ODE becomes :

with the sloe tends to -2/3

this can not solved analytically.

I am going to test the code for this equation now.

[sepideh]

xspan = [1 0];

y0 = [1e-8 -1/3];

fun = @(x,y)[y(2);-(3*y(2)^2+x*y(2)+y(1))/(3*y(1))];

[X,Y] = ode45(fun,xspan,y0);

plot(X,Y(:,1))

hold on

syms F(Y)

F(Y)= 1/6*(1-Y^2);

fplot(F,[0 1],'o')

axis([0 1 0 1])

hold on

xspan = [1 0];

y0 = [1e-8 -2/3];

fun = @(x,y)[y(2);-(3*y(2)^2+2*x*y(2)-y(1))/(3*y(1))];

[X,Y] = ode45(fun,xspan,y0);

plot(X,Y(:,1),'--')

hold off

as you can see the constant injection affect the height of shape function (match with our governing equation) , this is what I expected based on literature.

thank you

[sepideh]

Dear Torsten

I am trying to use the code for radial geometry but the numerical solution doesnot match the analytical one.

Would you please take a look at this and help me find out what is the problem? somehow in the mid interval its deviated

the ODE reads:

analytical solution is : (dots)

I used the base that you sent to me for backward method :

% plot analytical solution.

syms F(Y)

F(Y)= 1/8*(1-Y^2);

fplot(F,[0 1],'o')

axis([0 1 0 1])

hold on

% numerical solution for alpha=0

xspan = [1 0];

y0 = [1e-8 -1/4];

fun = @(x,y)[y(2);-(4*y(2)^2+(x)*y(2)+2*y(1))/(4*y(1))];

[X,Y] = ode45(fun,xspan,y0);

plot(X,Y(:,1))

ylabel('\Phi(\zeta)')

xlabel('\zeta')

hold off

[sepideh]

[Torsten]

syms x psi(x) f(x)
f = 1/8*(1-x^2);
ode = 0.5*psi+0.25*x*diff(psi,x)+diff(psi*diff(psi,x));
res = simplify(subs(ode,psi,f))

res(x) = 

[sepideh]

Thanks for your answer

I am trying to reproduce the results of an article, in that article the results are totally match.

sorry, I do not understand your solution, you want to show that there is a residual?

[Torsten]

Yes. The analytical expression 1/8*(1-x^2) is not a solution of your differential equation.

[sepideh]

yes you are right. I made a mistake in my modeling.

the correct ODE is:

now there is no residual.

but still there is a problem in the pragh.

xspan = [1 0];

y0 = [1e-10 -1/4];

fun = @(x,y)[y(2);-(2*x*y(1)+x^2*y(2)+4*x*y(2)^2+4*y(1)*y(2))/(4*x*y(1))];

[X,Y] = ode45(fun,xspan,y0);

plot(X,Y(:,1))

[Torsten]

Yo divide by x and y(1). Thus you have to start with something small and different from 0 for y(1) in the y0-vector and you must integrate up to something small and different from 0 for xspan(2) in the xspan-vector.

[sepideh]

thanks. now they are match