Coupeling of three equations BVP

Dear all,

I have a maths understanding problem with the following code. I tried to couple three ode equations with another and solved them with a bvp solver. As for the boundary conditions I applied all (just as an example) to the second derivative since in my problem I want to set Newmann boundary conditions.

However, what I noticed is that I am not able to give the third equation a dependence on C. The system is only solvable with Eq. dydx(6) decoupled from C.

How would I be able to solve the system with including a dependence of Eq. dydx(6) with C. Is a boundary condition to dydx(1), dydx(3) or dydx(5) nessecary in all cases if I want to solve this kind of system? I hope somebody more gifted in maths can help me with that.

Thanks for you help!

clear all;

close all;

clc;

xmesh = linspace(0, 1, 100);

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

sol = bvp4c(@ode, @bcfun, sol_init);

plot(sol.x, sol.y)

function dydx = ode(x,y)

A = y(1);

B = y(3);

C = y(5);

S = exp(A*B)-C;

dydx(1) = y(2);

dydx(2) = S;

dydx(3) = y(4);

dydx(4) = -S;

dydx(5) = y(6);

dydx(6) = exp(A+B);

% dydx(6) = C;

end

function res = bcfun(ya, yb);

res(1) = ya(2);

res(2) = yb(2);

res(3) = ya(4);

res(4) = yb(4);

res(5) = ya(6);

res(6) = yb(6);

end

Answer 1

[James Blanchard]

I get errors with both choices for dydx(6). They seem to result from your initial guesses (sol_init). Try something different.

[Thanh Hoang]

Hey sorry, yes you are right! I had a different expresion for dydx(6) originally:

dydx(6) = exp(A+B)

Here the corrected code

clear all;
close all;
clc;
xmesh = linspace(0, 1, 100);
sol_init = bvpinit(xmesh, [0 0 1 0 0 1]);
sol = bvp4c(@ode, @bcfun, sol_init);

Warning: Unable to meet the tolerance without using more than 1666 mesh points.
The last mesh of 892 points and the solution are available in the output argument.
The maximum residual is 0.0520983, while requested accuracy is 0.001.

plot(sol.x, sol.y)

function dydx = ode(x,y)
A = y(1);
B = y(3);
C = y(5);
S = exp(A*B)-C;
dydx(1) = y(2);
dydx(2) = S;
dydx(3) = y(4);
dydx(4) = -S;
dydx(5) = y(6);
dydx(6) = exp(A+B);
% dydx(6) = C;
end
function res = bcfun(ya, yb);
res(1) = ya(2);
res(2) = yb(2);
res(3) = ya(4);
res(4) = yb(4);
res(5) = ya(6);
res(6) = yb(6);
end

[Torsten]

There is no solution from bvp4c for

dydx(6) = exp(A+B);

either (see above).

The reason is that your equations only fix the derivative of y(5). So in my opinion, y(5) is only fixed by your equations up to an additive constant.