MATLAB: Solve a system of second order coupled differential equations with bvp4c

bvp4cdifferential equationsode45

Hi there, I am trying to use MatLabs bcp4c solver to solve the system:
eq 1: j'' = ( j' + t ) * (num1 – j + c'/c – t')
eq 2: c' = j' + t' – c
eq 3: t' = – j' * exp ( j)
with boundary conditions: j=1, c=0 at x=0 and j=0, c=1 at x=1.
I have linearized the problem, so that: Y(1) = j, Y(2)=j', Y(3) = c, Y(4) = t.
% Define space from 0 to 1.
x = linspace(0, 1, 1000);
% Inital guesses for the solver
guess = bvpinit(x, [1;0;1;1]);
% Solve function in space:
sol = bvp4c(@(x,Y) bvp_function(x,Y), @(Ya,Yb) boundary_conditons_fcn(Ya,Yb), guess);
% Export function:
solution_BVP_num=sol.y;
% Boundary conditions:
function res = boundary_conditons_fcn(Ya,Yb)
% j=1, c=0 at x=0 and j=0, c=1 at x=1.
    res = [Ya(1)-1; 
           Yb(1)-0
           Ya(3)-0
           Yb(3)-1]; 
end
function dYdx = bvp_function(x,Y)
% Y(1)=j, Y(2)=j', Y(3) = xsi, Y(4) = t 
    num1 = 2.331;
    dYdx(1) = Y(2);
    dYdx(2) = (Y(2) + Y(4)) * (num1 - Y(2) - dYdx(3)/Y(3) - dYdx(4));             
    dYdx(3) = Y(2) + dYdx(4) - Y(3);
    dYdx(4) = -Y(2)*exp(Y(1));
    dYdx = dYdx(:);
end
I get the error:
"Index exceeds matrix dimensions."
I think I have posed the question in a bad manner. Anybody knows how to solve this?
Hope so. 🙂

Best Answer

In the definition of bvp_function, when MATLAB reaches the line
dYdx(2) = (Y(2) + Y(4)) * (num1 - Y(2) - dYdx(3)/Y(3) - dYdx(4));
at that time, the values of dYdx(3) and dYdx(4) are undefined, so MATLAB throws an error.
You can solve this problem, either by changing the order of calculation like this
function dYdx = bvp_function(x, Y)
% Y(1)=j, Y(2)=j', Y(3) = xsi, Y(4) = t 

    num1 = 2.331;
    dYdx(4) = -Y(2)*exp(Y(1));
    dYdx(3) = Y(2) + dYdx(4) - Y(3);
    dYdx(1) = Y(2);
    dYdx(2) = (Y(2) + Y(4)) * (num1 - Y(2) - dYdx(3)/Y(3) + dYdx(4));
    dYdx = dYdx(:);
end
or replace the values of c' and t' in equation 1 and 2 like this
function dYdx = bvp_function(x, Y)
% Y(1)=j, Y(2)=j', Y(3) = xsi, Y(4) = t 
    num1 = 2.331;
    dYdx(1) = Y(2);
    dYdx(2) = (Y(2) + Y(4)) * (num1 - Y(2) - (Y(2) - Y(2)*exp(Y(1)) - Y(3))/Y(3) + Y(2)*exp(Y(1)));             
    dYdx(3) = Y(2) - Y(2)*exp(Y(1)) - Y(3);
    dYdx(4) = -Y(2)*exp(Y(1));
    dYdx = dYdx(:);
end
However, even after this, there seems to be a singularity in your system of ODEs. You need to analyze your equations to see what might be causing the singularity.
Related Question