Problems with a Multiple Boundary Value Problem

Hi everyone,

currently, I'm trying to find a solution for a Multiple Boundary Value Problem; in particular, I cannot find out a way to impose:

y'(0)=0 (where argument zero is the center of the entire interval, which is also the point xc inthe following code where I split in two separate intervals) :: the problem has been highlighted in the following code lines with "???" as a comment (bottom of the code, in bc function).

Here I attach the entire code:

clear all,clc,close

%% Define the mesh and initial guess

alpha=pi/6;

xc = 0; %define the midpoint where the additional b.c. has to be set

n_tot_nodes=100; %number of nodes to be used for overall computation

I=[-alpha alpha]; %overall interval where ODE solution is requested

xmesh_A=linspace(I(1),xc,floor(n_tot_nodes/2)); %first interval for ODE solution

xmesh_B=linspace(xc,I(2),floor(n_tot_nodes/2)); %second interval for ODE solution

xmesh=[xmesh_A xmesh_B]; %total interval for ODE solution (xc repreated)

yinit = [-0.2; 0.1; 0.3]; %define the initial vector

sol = bvpinit(xmesh,yinit); %define the initial guess on the mesh

sol = bvp5c(@(x,y,r) funzione(x,y,r), @bc, sol); %solve the BVP

%% LOCAL FUNCTIONS FOR MBVP SOLUTION

function dydx = funzione(x,y,region) % equations to be solved (in proper intervals)

dydx=zeros(3,1);

dydx(1)=y(1);

dydx(2)=y(2);

switch region %interval switch

case 1 % x in [0 xc]

dydx(3)=-(y(2).^5.*cos(x).^2.*-3.675715394304492e+17+y(2).^5.*sin(x).^2.*6.248716170317636e+18+y(2).*cos(x).^2.*(y(1).^2.*7.68e-4+y(2).^2./1.5625e+4).^(3.0./2.0).*4.34805803763114e+28-y(1).^2.*y(2).^3.*cos(x).^2.*8.82171694633078e+18-y(2).*sin(x).^2.*(y(1).^2.*7.68e-4+y(2).^2./1.5625e+4).^(3.0./2.0).*2.17402901881557e+29-y(1).^5.*cos(x).*sin(x).^2.*6.351636201358161e+20+y(1).^2.*y(2).^3.*sin(x).^2.*1.499691880876233e+20+y(2).^3.*y(3).^2.*sin(x).^2.*7.351430788608983e+17+y(2).^3.*cos(x).^2.*(y(1).^2.*7.68e-4+y(2).^2./1.5625e+4).^(3.0./2.0).*1.667553788025104e+19-y(2).^3.*sin(x).^2.*(y(1).^2.*7.68e-4+y(2).^2./1.5625e+4).^(3.0./2.0).*1.834309166827614e+20-y(1).^4.*y(2).*cos(x).^2.*5.293030167798468e+19+y(1).^5.*cos(x).*sin(x).*6.351636201358161e+20+y(1).^4.*y(2).*sin(x).^2.*1.534978748661556e+21+y(1).^2.*y(2).*cos(x).^2.*(y(1).^2.*7.68e-4+y(2).^2./1.5625e+4).^(3.0./2.0).*2.001064545630125e+20-y(1).^3.*cos(x).*sin(x).*(y(1).^2.*7.68e-4+y(2).^2./1.5625e+4).^(3.0./2.0).*3.601916182134225e+21+y(1).*cos(x).*sin(x).^2.*(y(1).^2.*7.68e-4+y(2).^2./1.5625e+4).^(3.0./2.0).*2.608834822578684e+29-y(1).^2.*y(2).*sin(x).^2.*(y(1).^2.*7.68e-4+y(2).^2./1.5625e+4).^(3.0./2.0).*1.420755827397389e+22-y(2).*y(3).^2.*sin(x).^2.*(y(1).^2.*7.68e-4+y(2).^2./1.5625e+4).^(3.0./2.0).*1.000532272815062e+20-y(1).^3.*y(2).^2.*cos(x).*sin(x).^2.*1.058606033559694e+20+y(1).*y(2).^4.*cos(x).*sin(x).*8.82171694633078e+18+y(1).^4.*y(3).*cos(x).*sin(x).*5.293030167798468e+19+y(2).^4.*y(3).*cos(x).*sin(x).*7.351430788608983e+17+y(1).*y(2).^3.*y(3).*sin(x).^2.*1.323257541949617e+19+y(1).^3.*y(2).*y(3).*sin(x).^2.*2.646515083899234e+20+y(1).^3.*cos(x).*sin(x).^2.*(y(1).^2.*7.68e-4+y(2).^2./1.5625e+4).^(3.0./2.0).*3.601916182134225e+21-y(1).*cos(x).*sin(x).*(y(1).^2.*7.68e-4+y(2).^2./1.5625e+4).^(3.0./2.0).*2.608834822578684e+29-y(3).*cos(x).*sin(x).*(y(1).^2.*7.68e-4+y(2).^2./1.5625e+4).^(3.0./2.0).*4.34805803763114e+28+y(1).^3.*y(2).^2.*cos(x).*sin(x).*1.58790905033954e+20-y(1).*y(2).^4.*cos(x).*sin(x).^2.*4.41085847316539e+18+y(1).^2.*y(2).*y(3).^2.*sin(x).^2.*1.323257541949617e+19-y(1).*y(2).*y(3).*sin(x).^2.*(y(1).^2.*7.68e-4+y(2).^2./1.5625e+4).^(3.0./2.0).*2.601383909319162e+21+y(1).^2.*y(2).^2.*y(3).*cos(x).*sin(x).*1.323257541949617e+19-y(1).*y(2).^2.*cos(x).*sin(x).*(y(1).^2.*7.68e-4+y(2).^2./1.5625e+4).^(3.0./2.0).*7.003725909705437e+20-y(1).^2.*y(3).*cos(x).*sin(x).*(y(1).^2.*7.68e-4+y(2).^2./1.5625e+4).^(3.0./2.0).*2.001064545630125e+20-y(2).^2.*y(3).*cos(x).*sin(x).*(y(1).^2.*7.68e-4+y(2).^2./1.5625e+4).^(3.0./2.0).*5.002661364075312e+19+y(1).*y(2).^2.*cos(x).*sin(x).^2.*(y(1).^2.*7.68e-4+y(2).^2./1.5625e+4).^(3.0./2.0).*3.001596818445187e+20)./(y(1).^4.*sin(x).^2.*5.293030167798468e+19+y(2).^4.*sin(x).^2.*7.351430788608983e+17-sin(x).^2.*(y(1).^2.*7.68e-4+y(2).^2./1.5625e+4).^(3.0./2.0).*4.34805803763114e+28+y(1).^2.*y(2).^2.*sin(x).^2.*1.323257541949617e+19-y(1).^2.*sin(x).^2.*(y(1).^2.*7.68e-4+y(2).^2./1.5625e+4).^(3.0./2.0).*2.001064545630125e+20-y(2).^2.*sin(x).^2.*(y(1).^2.*7.68e-4+y(2).^2./1.5625e+4).^(3.0./2.0).*5.002661364075312e+19);

case 2 % x in [xc 2]

dydx(3)=-(y(2).^5.*cos(x).^2.*-3.675715394304492e+17+y(2).^5.*sin(x).^2.*6.248716170317636e+18+y(2).*cos(x).^2.*(y(1).^2.*7.68e-4+y(2).^2./1.5625e+4).^(3.0./2.0).*4.34805803763114e+28-y(1).^2.*y(2).^3.*cos(x).^2.*8.82171694633078e+18-y(2).*sin(x).^2.*(y(1).^2.*7.68e-4+y(2).^2./1.5625e+4).^(3.0./2.0).*2.17402901881557e+29-y(1).^5.*cos(x).*sin(x).^2.*6.351636201358161e+20+y(1).^2.*y(2).^3.*sin(x).^2.*1.499691880876233e+20+y(2).^3.*y(3).^2.*sin(x).^2.*7.351430788608983e+17+y(2).^3.*cos(x).^2.*(y(1).^2.*7.68e-4+y(2).^2./1.5625e+4).^(3.0./2.0).*1.667553788025104e+19-y(2).^3.*sin(x).^2.*(y(1).^2.*7.68e-4+y(2).^2./1.5625e+4).^(3.0./2.0).*1.834309166827614e+20-y(1).^4.*y(2).*cos(x).^2.*5.293030167798468e+19+y(1).^5.*cos(x).*sin(x).*6.351636201358161e+20+y(1).^4.*y(2).*sin(x).^2.*1.534978748661556e+21+y(1).^2.*y(2).*cos(x).^2.*(y(1).^2.*7.68e-4+y(2).^2./1.5625e+4).^(3.0./2.0).*2.001064545630125e+20-y(1).^3.*cos(x).*sin(x).*(y(1).^2.*7.68e-4+y(2).^2./1.5625e+4).^(3.0./2.0).*3.601916182134225e+21+y(1).*cos(x).*sin(x).^2.*(y(1).^2.*7.68e-4+y(2).^2./1.5625e+4).^(3.0./2.0).*2.608834822578684e+29-y(1).^2.*y(2).*sin(x).^2.*(y(1).^2.*7.68e-4+y(2).^2./1.5625e+4).^(3.0./2.0).*1.420755827397389e+22-y(2).*y(3).^2.*sin(x).^2.*(y(1).^2.*7.68e-4+y(2).^2./1.5625e+4).^(3.0./2.0).*1.000532272815062e+20-y(1).^3.*y(2).^2.*cos(x).*sin(x).^2.*1.058606033559694e+20+y(1).*y(2).^4.*cos(x).*sin(x).*8.82171694633078e+18+y(1).^4.*y(3).*cos(x).*sin(x).*5.293030167798468e+19+y(2).^4.*y(3).*cos(x).*sin(x).*7.351430788608983e+17+y(1).*y(2).^3.*y(3).*sin(x).^2.*1.323257541949617e+19+y(1).^3.*y(2).*y(3).*sin(x).^2.*2.646515083899234e+20+y(1).^3.*cos(x).*sin(x).^2.*(y(1).^2.*7.68e-4+y(2).^2./1.5625e+4).^(3.0./2.0).*3.601916182134225e+21-y(1).*cos(x).*sin(x).*(y(1).^2.*7.68e-4+y(2).^2./1.5625e+4).^(3.0./2.0).*2.608834822578684e+29-y(3).*cos(x).*sin(x).*(y(1).^2.*7.68e-4+y(2).^2./1.5625e+4).^(3.0./2.0).*4.34805803763114e+28+y(1).^3.*y(2).^2.*cos(x).*sin(x).*1.58790905033954e+20-y(1).*y(2).^4.*cos(x).*sin(x).^2.*4.41085847316539e+18+y(1).^2.*y(2).*y(3).^2.*sin(x).^2.*1.323257541949617e+19-y(1).*y(2).*y(3).*sin(x).^2.*(y(1).^2.*7.68e-4+y(2).^2./1.5625e+4).^(3.0./2.0).*2.601383909319162e+21+y(1).^2.*y(2).^2.*y(3).*cos(x).*sin(x).*1.323257541949617e+19-y(1).*y(2).^2.*cos(x).*sin(x).*(y(1).^2.*7.68e-4+y(2).^2./1.5625e+4).^(3.0./2.0).*7.003725909705437e+20-y(1).^2.*y(3).*cos(x).*sin(x).*(y(1).^2.*7.68e-4+y(2).^2./1.5625e+4).^(3.0./2.0).*2.001064545630125e+20-y(2).^2.*y(3).*cos(x).*sin(x).*(y(1).^2.*7.68e-4+y(2).^2./1.5625e+4).^(3.0./2.0).*5.002661364075312e+19+y(1).*y(2).^2.*cos(x).*sin(x).^2.*(y(1).^2.*7.68e-4+y(2).^2./1.5625e+4).^(3.0./2.0).*3.001596818445187e+20)./(y(1).^4.*sin(x).^2.*5.293030167798468e+19+y(2).^4.*sin(x).^2.*7.351430788608983e+17-sin(x).^2.*(y(1).^2.*7.68e-4+y(2).^2./1.5625e+4).^(3.0./2.0).*4.34805803763114e+28+y(1).^2.*y(2).^2.*sin(x).^2.*1.323257541949617e+19-y(1).^2.*sin(x).^2.*(y(1).^2.*7.68e-4+y(2).^2./1.5625e+4).^(3.0./2.0).*2.001064545630125e+20-y(2).^2.*sin(x).^2.*(y(1).^2.*7.68e-4+y(2).^2./1.5625e+4).^(3.0./2.0).*5.002661364075312e+19);

end

end

%-------------------------------------------

function res = bc(YL,YR) % boundary conditions (in proper intervals)

Q=1e-9;

res = [YL(1,1) % y_1_sx_(-alpha)= 0

------- %???

YR(1,1) - YL(1,2) % Continuity of y_1_(x) at x=0

YR(2,1) - YL(2,2) % Continuity of y_2_(x) at x=0

YR(3,1) - YL(3,2) % Continuity of y_3_(x) at x=0

YR(1,2)]; % y_1_dx_(alpha) = 0

end

%-------------------------------------------

Thanks in advance to everyone that will help with this problem. Best regards.

Answer 1

[Abhinaya Kennedy]

Hi Alessio,

To impose the condition (y'(0) = 0) for a boundary value problem in MATLAB, you can modify the boundary condition function “bc” to include this derivative condition. Assuming “y'(0) = 0” applies to the first component of your solution vector “y” and considering you're dealing with continuity at “xc = 0”, here's an altered version of your “bc” function:

function res = bc(YL,YR)
    res = [YL(1);             % Left boundary condition for y    
        YL(2);             % y'(0) = 0, assuming YL(2) represents this derivative
        YR(1) - YL(1);     % Continuity of y at x=0
        YR(2) - YL(2);     % Continuity of y' at x=0, if applicable
        YR(3) - YL(3);     % Additional continuity conditions, if any
        YR(1)];            % Right boundary condition for y
end

This assumes:

• “YL(2)” represents (y'(0)) from the left side. Adjust according to your problem's specifics and how your differential equations are structured.
• Continuity conditions and boundary conditions are correctly defined for your problem's context.

This link will give you more insights into solving boundary value problems with multiple boundary conditions:

https://www.mathworks.com/help/matlab/math/solve-bvp-with-multiple-boundary-conditions.html

Hope this helps!