How to use correctly bvp4c / bvp5c

Hello,

How can I solve this boundaries value problem : y''(r) = (k*(k+1)/r^2)*y(r)-((2*k+1)/r)*f_1(r)*f_2(r)

where :

1. r is defined as r=0:dr:r_max with dr and r_max given.
2. y(r) is the function that I need to solve.
3. k is a given integer parameter.
4. f_1(r) and f_2(r) are two known functions.

The boundaries conditions are the following:

1. y(0) = 0
2. y'(r_max) = -(k/r_max)*y(r_max)

I have tried solving it using bvp5c but I failed :

 if
function sol=myproblem(dr,r_max,k,f1,f2)
r=0:dr:r_max;
solinit=bvpinit([0 r_max],[0 0]);
sol=bvp5c(@odefun,@bcfun,solinit);
end

 if
function dydr=odefun(r,y,k,f1,f2,dr)
c1=k*(k+1)./r.^2;
c1(1)=k*(k+1)./(dr/10)^2; % to avoid dividing by 0
c2=(2*k+1)*f1(:).*f2(:)./r;
c2(1)=0; % to avoid dividing by 0
dydr=[y(2);c1.*y(1)-c2];
end

if
function res=bcfun(ya,yb,k,r_max)
res=[ya(1);yb(2)+(k/r_max)*yb(1)];
end

In the code above, ignore the if, I can not remove them without loosing the code format.

thanks in advance,

David

Answer 1

[David Liegeois]

Hi Torsten,

Thank you for your answer.

Somme addings regarding what I wrote : f1 and f2 have already the same length as r when I pass them into the function myproblem. I do not need to make an interpolation.

Matlab writes : Error in bvp5c : Not enough input argument.

David

[Torsten]

Please show the code you are using at the moment.

Best wishes

Torsten.

[David Liegeois]

 1
z=1;
n=[1 2];
dr=0.01;
r_max=10;
k=1;
r=0:dr:r_max;
x(1,:)=z*r./n(1);
x(2,:)=z*r./n(2);
p1=exp(-x(1,:)).*x(1,:).^n(1);
p2=exp(-x(2,:)).*x(2,:).^n(2);
clear x z n;
sol=myproblem(dr,r_max,k,p1,p2);

 2
function sol=myproblem(dr,r_max,k,p1,p2)
r=0:dr:r_max;
c1=k.*(k+1)./r.^2;
c1(1)=k.*(k+1)/(dr/10)^2; % to avoid dividing by 0
c2=(2*k+1)*p1.*p2./r;
c2(1)=0; % to avoid dividing by 0
c=k/r_max;
solinit=bvpinit(r,p1.*p2);
sol=bvp5c(@(r,y)odefun(c1,c2),@(ya,yb)bcfun(c),solinit);
end

 3
function dydr=odefun(r,y,c1,c2)
dydr=[y(2);c1.*y(1)-c2];
end

 4
function res=bcfun(ya,yb,c)
res=[ya(1);yb(2)+c*yb(1)];
end

[Torsten]

Try

1
z=1;
n=[1 2];
dr=0.01;
r_max=10;
k=1;
r=0:dr:r_max;
x(1,:)=z*r./n(1);
x(2,:)=z*r./n(2);
p1=exp(-x(1,:)).*x(1,:).^n(1);
p2=exp(-x(2,:)).*x(2,:).^n(2);
clear x z n;
sol=myproblem(dr,r_max,k,p1,p2);

2
function sol=myproblem(dr,r_max,k,p1,p2)
r=0:dr:r_max;
c1=k.*(k+1)./r.^2;
c1(1)=k.*(k+1)/(dr/10)^2; % to avoid dividing by 0
c2=(2*k+1)*p1.*p2./r;
c2(1)=0; % to avoid dividing by 0
c=k/r_max;
solinit=bvpinit(r,[1 1]);
sol = bvp5c(@(r,y)odefun(r,y,r,c1,c2),@(ya,yb)bcfun(ya,yb,c),solinit)

3
function dydr = odefun(r,y,rs,c1,c2)
c1_actual=interp1(rs,c1,r);
c2_actual=interp1(rs,c2,r);
dydr = [y(2);c1_actual*y(1)-c2_actual];
end

4
function res=bcfun(ya,yb,c)
res=[ya(1);yb(2)+c*yb(1)];
end

Best wishes

Torsten.

[David Liegeois]

Hi Torsten,

I try it but it fails with the error : The derivative function ODEFUN should return a column vector of length 2.

I think the problem comes from the structure solinit. Inside it, I have :

1. solver : 'bvpinit'
2. x : 1x1001 double
3. y : 2x1001 double
4. yinit : [0 0]

Following the error message, x and y should be 1001x2 double but I can not succeed to transpose them.

David

[Torsten]

Try

2
function sol=myproblem(dr,r_max,k,p1,p2)
r=0:dr:r_max;
c1=k.*(k+1)./r.^2;
c1(1)=k.*(k+1)/(dr/10)^2; % to avoid dividing by 0
c2=(2*k+1)*p1.*p2./r;
c2(1)=0; % to avoid dividing by 0
c=k/r_max;
solinit=bvpinit(r,[1 1]);
sol = bvp5c(@(rr,y)odefun(rr,y,r,c1,c2),@(ya,yb)bcfun(ya,yb,c),solinit)

Best wishes

Torsten.

[David Liegeois]

Hi Torsten.

Thank you, t works (I have an output) but the problem is that I gave a wrong solinit. With your solution I have a yinit constant (=1 for y(1,:) and y(2,:)) but I would like to give a function depending of r. When I try this :

 1
 ...
yinit(1,:)=p1.*p2;
yinit(2,:)=-k*yinit(1,:)./r;
yinit(2,1)=0;
solinit=bvpinit(r,yinit);
...

I have an error : The guess for the solution must return a vector.

Do you have an idea ?

David

[Torsten]

Take a look at the example

"Compute Fourth Eigenvalue of Mathieu’s Equation"

under

https://de.mathworks.com/help/matlab/ref/bvp4c.html

Best wishes

Torsten.

[David Liegeois]

Hello Torsten,

I have tried to reproduce the same in my code but I failed. With 1 in function my problem before the call to bvp5c and with the function mat4init, I have the error : Not enough input arguments.

 1
c3=r.*P(1,:).*P(1,:);
c4=-k*c3./r;
c4(1)=0;
solinit=bvpinit(r,@mat4init);

 2
function yinit=mat4init(c3,c4)
yinit=[c3;c4];
end

David

[Torsten]

1
c3=r.*P(1,:).*P(1,:);
c4=-k*c3./r;
c4(1)=0;
solinit=bvpinit(r,@(rr)mat4init(rr,r,c3,c4));

2
function yinit=mat4init(rr,r,c3,c4)
c3_actual=interp1(r,c3,rr);
c4_actual=interp1(r,c4,rr);
yinit=[c3_actual;c4_actual];
end

Best wishes

Torsten.

[David Liegeois]

Ok everything works ! Thanks again Torsten.

Answer 2

[Torsten]

1.

Build vectors r1 and r2 in the calling program which are of the same length as f1 and f2 and which contain the values of r at which f1 and f2 are given.

r1=0:dr:r_max 
r2=0:dr:r_max

2.

Call bvp5c as

sol = bvp5c(@(r,y)odefun(r,y,k,r1,f1,r2,f2),@(ya,yb)bcfun(ya,yb,k,r_max),solinit)

3.

Use the following function as "odefun":

function dydr = odefun(r,y,k,r1,f1,r2,f2)
c1 = k*(k+1)/r^2;  
f1_actual=interp1(r1,f1,r);
f2_actual=interp1(r2,f2,r);
c2 = (2*k+1)/r*f1_actual*f2_actual;  
dydr = [y(2);c1*y(1)-c2];
end

Best wishes

Torsten.

[Colin Wang]

Hi Torsen,

Thank you for your answer! Use this example I'm able to solve my equations and get correct answer.

However there is one part I dont understand, really appreciate it if you can teach me by explain why it doesnt work if written as:

 sol = bvp5c(@(r,y)odefun(r,y,r,c1,c2),@(ya,yb)bcfun(ya,yb,c),solinit)

and becomes working written as :

 sol = bvp5c(@(rr,y)odefun(rr,y,r,c1,c2),@(ya,yb)bcfun(ya,yb,c),solinit)

'rr' is not defined in the code, and in your code:

 3
function dydr = odefun(r,y,rs,c1,c2)
c1_actual=interp1(rs,c1,r);
c2_actual=interp1(rs,c2,r);
dydr = [y(2);c1_actual*y(1)-c2_actual];
end

the parameters in odefun is r,y,rs,c1,c2

but called as odefun(rr,y,r,c1,c2). where is rs? and why can the parameters changes from (r,y,rs,c1,c2) to (rr,y,r,c1,c2)?

Thank you in advance for your helpful answers!!

Rgds,

Colin

[Colin Wang]

Hi Torsen,

I tried to undertand it, here are some thoughts:

 sol = bvp5c(@(r,y)odefun(r,y,r,c1,c2),@(ya,yb)bcfun(ya,yb,c),solinit)

in this sentence, it doesnt work bcuz the 1st and 3rd parameter in odefun are both 'r', which is duplicated and wrong, as the 1st 'r' is not involved in the calculation, it can be written as 'rr' or anything just to avoid the duplication to the 3rd parameter 'r'.

and for the inconsistency define of parameters bewteen the 1st parameter ('rr') at calling the bvp5c

 sol = bvp5c(@(rr,y)odefun(rr,y,r,c1,c2),@(ya,yb)bcfun(ya,yb,c),solinit)

and part 3, the 1st parameter ('r') in the definition of odefun

 3
function dydr = odefun(r,y,rs,c1,c2)
c1_actual=interp1(rs,c1,r);
c2_actual=interp1(rs,c2,r);
dydr = [y(2);c1_actual*y(1)-c2_actual];
end

Same reason, as the 1st paramater is not involved in the calculation, so the incosistency for the 1st parameter doesnt matter?

However I still dont undertand the rs thing. The interp1 should do its work on (rs, c1/c2), but rs is not defined?

As provided, f1 and f2 have already the same length as r , which means rs=r? but there is no definition of this, how can the code gives the result?

I think it can only works if use rs replaces all the r in:

2
function sol=myproblem(dr,r_max,k,p1,p2)
r=0:dr:r_max;
c1=k.*(k+1)./r.^2;
c1(1)=k.*(k+1)/(dr/10)^2; % to avoid dividing by 0
c2=(2*k+1)*p1.*p2./r;
c2(1)=0; % to avoid dividing by 0
c=k/r_max;
solinit=bvpinit(r,[1 1]);
sol = bvp5c(@(rr,y)odefun(rr,y,r,c1,c2),@(ya,yb)bcfun(ya,yb,c),solinit)

or predefine rs=r.

Kindly let me know ur comments.

Rgds,

Colin

[Colin Wang]

Hi Torsen,

I've gone throuth the page of 'Parameterizing Functions'

https://ww2.mathworks.cn/help/matlab/math/parameterizing-functions.html?lang=en

I think the answer of my questions is related with this?

Although not mentioned in above page, I tried with a simple function:

b = 2;
c = 3.5;
x = fzero(@(x) cubicpoly(x,b,c),0)
function y = cubicpoly(x,b,c)
y = x^3 + b*x + c;
end

and even I changed the code to

d = 2;
f = 3.5;
x = fzero(@(xxx) cubicpoly(xxx,d,f),0)
function y = cubicpoly(x,b,c)
y = x^3 + b*x + c;
end

It gives the same result. This verified my guess, that in your code, u wrote

odefun(r,y,rs,c1,c2)

in part 3 definition but wrote as

sol = bvp5c(@(rr,y)odefun(rr,y,r,c1,c2),@(ya,yb)bcfun(ya,yb,c),solinit)

to call the odefun. Here 'rs' has the same value with 'r' and it doesnt matter if you write the 1st parameter as 'rr' or 'r', they are just symbols to present the parameter?

However I couldnt find a document from Matlab describing more details on it, only the above link.

Rgds,

Colin