Syntax of res-function for bvp4c

function res = bcfun(ya,yb,I)
ya(1)=850;
yb(1)=303;
ya(3)=850;
yb(3)=303;
ya(5)=r_cu;      %constant
yb(5)=r_cu-((rho_p_T_c/A_p+rho_n_T_c/A_n)-(rho_p_T_h/A_p+rho_n_T_h/A_n))*L;            %all constants
V_tho=2*I*(r_cu+(rho_p_T_c/A_p+rho_n_T_c/A_n)*L)-(alpha_p_T_h-alpha_n_T_h)*T_h-(alpha_p_T_c-alpha_n_T_c)*T_c; %all constants except for parameter I

res = [ya(1)-yb(1)
         ya(3)-yb(1)
             ya(2)
             ya(4)
             ya(5)
             yb(5)
             ya(6)
            yb(6)-V_tho];
end

[Torsten]

If you want to set yi at x=0 to a value "v", res = ya(i) - v.

If you want to set yi at x=L to a value "v", res = yb(i) - v.

[Justin Riggio]

Thanks Torsten for the reponse. I am still a bit stuck. Using your suggestion, there would be 10 boundary conditions when the formulation has 6 variables and 1 parameter. Do you think there is any other way to prescribe this many?

[Torsten]

If you have six ODEs and one parameter, you need seven boundary conditions. And bvp4c will complain if you prescribe more or less.

[Justin Riggio]

Yes. I have tried reducing the number of boundary conditions and depending on which ones used, the solution diverges or has an error because of a singular Jacobian. I will have to satisfy bvp4c someway or another and looks like it has to be by using more ODEs. Another question I have that can help solve this has to do with the possibility of having a boundary condition that is a function of a variable at one position. In this case, both ya(2) and ya(4), is a function of y(1) at x=1 and y(3) at x=1 respectively. Have you ever delt with this?

[Torsten]

So y(2) at x=0 is a function of y(1) at x=1 and y(4) at x=0 is a function of y(3) at x=1 ? This is only possible if an explicit value for y(2) at x=0 or y(1) at x=1 is prescribed as a second boundary condition for the first case and if an explicit value for y(4) at x=0 or y(3) at x=1 is prescribed as a second boundary condition for the second case.

[Justin Riggio]

Yes exactly.

This is what I should have:

ya(2)+k_p_T_h*A_p*(T_h-y(1)); %how can I tell it to use y(1) at x=1
ya(4)+k_n_T_h*A_p*(T_h-y(3)); %how can I tell it to use y(3) at x=1

With this, I am not sure how to prescribe another boundary condition for both y(2) at x=0 and y(4) at x=0.

Perhaps I can only prescribe y(1) at x=1 and y(3) at x=1 (have yet to look into how to do so), however they need to be changing as the solution converges.

[Torsten]

If x=1 is the right boundary point, you can set

res = ya(2)-k_p_t_h*A_p*(T_h-yb(1))

But ya(2) or yb(1) must be specified somewhere else as explicit boundary condition, e.g.

res = ya(2) - 20

If x=1 is not a boundary point, the condition cannot be set in bvp4c.

[Justin Riggio]

x=1 is not a boundary point. There is more to work out than just syntax. Thanks Torsten

[Justin Riggio]

Hi Torsten,

In the res function, I would like to put guess values at specific points into those boundary conditions as such:

[ya(2)-yb(2)-(k_p_T_h*A_p*(T_h-solinit.y(1,2)/solinit.x(1,2))-k_p_T_c*A_p*(solinit.y(1,9)-T_c)/...(solinit.x(10)-solinit.x(9)))]

I am getting an error: undefined variable or class "solinit.y"

Is there a way this can work? I am thinking of using continuation afterwards

Justin

[Torsten]

Call bvp4c like

sol = bvp4c(@odefun,@(ya,yb)bcfun(ya,yb,solinit),solinit);

to make "solinit" available in "bcfun".

[Justin Riggio]

I just tried and I got the error:

"Error using main>@(ya,yb)bcfun(ya,yb,solinit)

Too many input arguments."

[Torsten]

Did you also change the argument list of the bcfun-routine itself:

function res = bcfun(ya,yb,solinit)

?

[Justin Riggio]

Yes,

function res = bcfun(ya,yb,I,solinit)

I is a parameter

[Torsten]

Then

sol = bvp4c(@odefun,@(ya,yb,l)bcfun(ya,yb,l,solinit),solinit);

[Justin Riggio]

I did so and have yet again too many input arguments.

This is how the res function looks now:

function res = bcfun(ya,yb,I,solinit)
constants
res = [ya(1)-yb(1)-547;...
       ya(2)-yb(2)-(k_p_T_h*A_p*(T_h-solinit.y(1,2)/solinit.x(1,2))-k_p_T_c*A_p*(solinit.y(1,9)-T_c)/(solinit.x(1,10)-solinit.x(1,9)));...
       ya(3)-yb(3)-547;...
       ya(4)-yb(4)-(k_n_T_h*A_n*(T_h-solinit.y(3,2)/solinit.x(1,2))-k_n_T_c*A_n*(solinit.y(3,9)-T_c)/(solinit.x(1,10)-solinit.x(1,9)));...
       ya(5)-yb(5)-((rho_p_T_c/A_p+rho_n_T_c/A_n)-(rho_p_T_h/A_p+rho_n_T_h/A_n))*L;...
       ya(6)
       yb(6)-V_tho];
end 

[Torsten]

Try if it works if you use the line

global solinit

in the program where you call bvp4c as well as in bcfun.

[Justin Riggio]

Seems to work! Although I have gotten the ''unable to solve the collocation equations -- a singular Jacobian encountered.'' error, there was an attempt at least. I will try to implement continuation, however not sure that it will help this error

[Torsten]

I think the error before arises because you did not initialize your parameter l within the solinit construct.

Take a look at the example

"Compute Fourth Eigenvalue of Mathieu’s Equation"

under

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

especially the definition of "solinit".

[Justin Riggio]

I have been doing so this way:

solinit=bvpinit(x,@y_init,0.001);

[Justin Riggio]

Also when I implement the global solinit I had to change the sol construct back to the original:

sol=bvp4c(@odefun,@bcfun,solinit);
As there were too many input arguments with
sol = bvp4c(@odefun,@(ya,yb,l)bcfun(ya,yb,l,solinit),solinit);

[Justin Riggio]

Hi Torsten, I am attempting to use a while loop for continuation like this which is obviously incorrect. I am wondering of a way I can store sol.y from the solver and compare it to the previous value for this case. Perhaps there is a better approach.

sol=bvpinit(x,@y_init,I);
while abs(sol.y(1,2) - sol.y(1,2)) > 0.000001  
sol=bvp4c(@odefun,@bcfun,sol);
end

[madhan ravi]

@Justin i edited all your comment next time format it properly so that it's easy to read

[Justin Riggio]

Thank you!

[Justin Riggio]

sol=bvpinit(x,@y_init,I);
old_value=5;

while abs(sol.y(1,2) - old_value) > 0.000001 

sol=bvpinit(x,@y_init,I);
sol=bvp5c(@odefun,@bcfun,sol);
sol.y(1,2)=old_value;
end

[Justin Riggio]

I have gotten to this above, however, when I check the res function, it does not change then reaches the error once again for the encounter of a singular Jacobian!