Error solving bvp4c - Singular jacobian

Hello everyone.

I'm trying to solve in Matlab2017b an ODE with the boundary conditions:

,

For this purpose, I have used the solver bvp4c. I think that this equation must be solvable for all values of [z1 z2 z3] because it has the form of a generic forced oscillator. However, there are many for which appears the error: singular jacobian (like the values I have written down) and I cannot guess which is the problem. Any idea?

Thank you in advance!

%Constants
hb = 6.626e-34/(2*pi);
m = 9*1.660538921e-27;
w0 = 2e6*2*pi;
Cc = (1.6e-19)^2/(4*pi*8.854e-12);
R = 10;
gammam = 1;
gammap = (3*R^3/(R^3+2))^(1/4);
NT = 50;
n = 100;
tf = 3.2e-6;
%Initial seed
z=[-104.2545  628.8529   33.2914];
z1=z(1); z2 =z(2); z3=z(3);

New1 = bvp4c(@(t,y) new_qubic(t, y, z1, z2, z3, gammam, gammap, tf, w0, Cc, m),@bvp4bc,solinit,options);

function dydx=new_qubic(t, y, z1, z2, z3, gammam, gammap, tf, w0, Cc, m)

z4=0;
    
    rhop=1+(-126-z1-5*z2-15*z3-35*z4+126*gammap)*t^5/tf^5+...
        (420+5*z1+24*z2+70*z3+160*z4-420*gammap)*t^6/tf^6+...
        (-540-10*z1-45*z2-126*z3-280*z4+540*gammap)*t^7/tf^7+...
        (315+10*z1+40*z2+105*z3+224*z4-315*gammap)*t^8/tf^8+...
        (-70-5*z1-15*z2-35*z3-70*z4+70*gammap)*t^9/tf^9+...
        z1*t^10/tf^10+z2*t^11/tf^11+z3*t^12/tf^12+z4*t^13/tf^13;
    rho1p=5*(-126-z1-5*z2-15*z3-35*z4+126*gammap)*t^4/tf^5+...
         6*(420+5*z1+24*z2+70*z3+160*z4-420*gammap)*t^5/tf^6+...
         7*(-540-10*z1-45*z2-126*z3-280*z4+540*gammap)*t^6/tf^7+...
         8*(315+10*z1+40*z2+105*z3+224*z4-315*gammap)*t^7/tf^8+...
         9*(-70-5*z1-15*z2-35*z3-70*z4+70*gammap)*t^8/tf^9+...
         10*z1*t^9/tf^10+11*z2*t^10/tf^11+12*z3*t^11/tf^12+13*z4*t^12/tf^13;
    rho2p=20*(-126-z1-5*z2-15*z3-35*z4+126*gammap)*t^3/tf^5+...
         30*(420+5*z1+24*z2+70*z3+160*z4-420*gammap)*t^4/tf^6+...
         42*(-540-10*z1-45*z2-126*z3-280*z4+540*gammap)*t^5/tf^7+...
         56*(315+10*z1+40*z2+105*z3+224*z4-315*gammap)*t^6/tf^8+...
         72*(-70-5*z1-15*z2-35*z3-70*z4+70*gammap)*t^7/tf^9+...
         90*z1*t^8/tf^10+110*z2*t^9/tf^11+132*z3*t^10/tf^12+156*z4*t^11/tf^13;
    rho3p=60*(-126-z1-5*z2-15*z3-35*z4+126*gammap)*t^2/tf^5+...
         120*(420+5*z1+24*z2+70*z3+160*z4-420*gammap)*t^3/tf^6+...
         210*(-540-10*z1-45*z2-126*z3-280*z4+540*gammap)*t^4/tf^7+...
         336*(315+10*z1+40*z2+105*z3+224*z4-315*gammap)*t^5/tf^8+...
         504*(-70-5*z1-15*z2-35*z3-70*z4+70*gammap)*t^6/tf^9+...
         720*z1*t^7/tf^10+990*z2*t^8/tf^11+1320*z3*t^9/tf^12+1716*z4*t^10/tf^13;
    rho4p=120*(-126-z1-5*z2-15*z3-35*z4+126*gammap)*t^1/tf^5+...
         360*(420+5*z1+24*z2+70*z3+160*z4-420*gammap)*t^2/tf^6+...
         840*(-540-10*z1-45*z2-126*z3-280*z4+540*gammap)*t^3/tf^7+...
         1680*(315+10*z1+40*z2+105*z3+224*z4-315*gammap)*t^4/tf^8+...
         3024*(-70-5*z1-15*z2-35*z3-70*z4+70*gammap)*t^5/tf^9+...
         5040*z1*t^6/tf^10+7920*z2*t^7/tf^11+11880*z3*t^8/tf^12+17160*z4*t^9/tf^13;

    rhom=126*(gammam-1)*t^5/tf^5-420*(gammam-1)*t^6/tf^6+...
          540*(gammam-1)*t^7/tf^7-315*(gammam-1)*t^8/tf^8+70*(gammam-1)*t^9/tf^9+1;
    rho1m=630*(gammam-1)*t^4/tf^5-2520*(gammam-1)*t^5/tf^6+...
           3780*(gammam-1)*t^6/tf^7-2520*(gammam-1)*t^7/tf^8+630*(gammam-1)*t^8/tf^9;
    rho2m=2520*(gammam-1)*t^3/tf^5-12600*(gammam-1)*t^4/tf^6+...
           22680*(gammam-1)*t^5/tf^7-17640*(gammam-1)*t^6/tf^8+5040*(gammam-1)*t^7/tf^9;
    rho3m=7560*(gammam-1)*t^2/tf^5-50400*(gammam-1)*t^3/tf^6+...
           113400*(gammam-1)*t^4/tf^7-105840*(gammam-1)*t^5/tf^8+35280*(gammam-1)*t^6/tf^9;
    rho4m=15120*(gammam-1)*t^1/tf^5-151200*(gammam-1)*t^2/tf^6+...
           453600*(gammam-1)*t^3/tf^7-529200*(gammam-1)*t^4/tf^8+211680*(gammam-1)*t^5/tf^9;

    wp=sqrt((sqrt(3)*w0)^2/rhop^4-rho2p/rhop);
    w1p=1/2/wp*(-4*(sqrt(3)*w0)^2*rho1p/rhop^5-(rho3p*rhop-rho2p*rho1p)/rhop^2);
    w2p=-1/4/wp^3*(-4*(sqrt(3)*w0)^2*rho1p/rhop^5-(rho3p*rhop-rho2p*rho1p)/rhop^2)^2+...
        1/(2*wp)*(-(4*(sqrt(3)*w0)^2*rho2p*rhop-20*(sqrt(3)*w0)^2*rho1p^2)/rhop^6-(rho4p*rhop^2-rho2p^2*rhop-2*rho3p*rho1p*rhop-2*rho2p*rho1p^2)/rhop^3);

    wm=sqrt(w0^2/rhom^4-rho2m/rhom);
    w1m=1/2/wm*(-4*w0^2*rho1m/rhom^5-(rho3m*rhom-rho2m*rho1m)/rhom^2);
    w2m=-1/4/wm^3*(-4*w0^2*rho1m/rhom^5-(rho3m*rhom-rho2m*rho1m)/rhom^2)^2+...
        1/(2*wm)*(-(4*w0^2*rho2m*rhom-20*w0^2*rho1m^2)/rhom^6-(rho4m*rhom^2-rho2m^2*rhom-2*rho3m*rho1m*rhom-2*rho2m*rho1m^2)/rhom^3);

    ddd=(4*2^(2/3)*Cc^(1/3)*(-2*m*wm*w1m+2*m*wp*w1p)^2)/(9*(-m*wm^2+m*wp^2)^(7/3))-...
        (2^(2/3)*Cc^(1/3)*(-2*m*w1m^2+2*m*w1p^2-2*m*wm*w2m+2*m*wp*w2p))/(3*(-m*wm^2+m*wp^2)^(4/3));

    dydx=[y(2) -sqrt(m/2)*ddd-wp^2*y(1)];

end

function res = bvp4bc(ya,yb)
   res = [ya(1) ya(2)];
end

[Torsten]

Neither "solinit" nor "options" is supplied.

[Jesús Parejo]

ups, sorry, these are solinit and options

nt= 2000;
tf=3.2e-6;
solinit = bvpinit(linspace(0,tf,nt),[0 0]);
options = bvpset('RelTol',10^(-6));

Answer 1

[Torsten]

Try this code following John's suggestion:

%Constants
hb = 6.626e-34/(2*pi);
m = 9*1.660538921e-27;
w0 = 2e6*2*pi;
Cc = (1.6e-19)^2/(4*pi*8.854e-12);
R = 10;
gammam = 1;
gammap = (3*R^3/(R^3+2))^(1/4);
NT = 50;
n = 100;
tf = 3.2e-6;
%Initial seed
z=[-104.2545  628.8529   33.2914];
z1=z(1); z2 =z(2); z3=z(3);
[T,Y]=ode15s(@(t,y) new_qubic(t, y, z1, z2, z3, gammam, gammap, tf, w0, Cc, m),[0 tf],[0 0]);
plot(T,Y(:,1))

function dydx=new_qubic(t, y, z1, z2, z3, gammam, gammap, tf, w0, Cc, m)
z4=0;

rhop=1+(-126-z1-5*z2-15*z3-35*z4+126*gammap)*t^5/tf^5+...
    (420+5*z1+24*z2+70*z3+160*z4-420*gammap)*t^6/tf^6+...
    (-540-10*z1-45*z2-126*z3-280*z4+540*gammap)*t^7/tf^7+...
    (315+10*z1+40*z2+105*z3+224*z4-315*gammap)*t^8/tf^8+...
    (-70-5*z1-15*z2-35*z3-70*z4+70*gammap)*t^9/tf^9+...
    z1*t^10/tf^10+z2*t^11/tf^11+z3*t^12/tf^12+z4*t^13/tf^13;
rho1p=5*(-126-z1-5*z2-15*z3-35*z4+126*gammap)*t^4/tf^5+...
    6*(420+5*z1+24*z2+70*z3+160*z4-420*gammap)*t^5/tf^6+...
    7*(-540-10*z1-45*z2-126*z3-280*z4+540*gammap)*t^6/tf^7+...
    8*(315+10*z1+40*z2+105*z3+224*z4-315*gammap)*t^7/tf^8+...
    9*(-70-5*z1-15*z2-35*z3-70*z4+70*gammap)*t^8/tf^9+...
    10*z1*t^9/tf^10+11*z2*t^10/tf^11+12*z3*t^11/tf^12+13*z4*t^12/tf^13;
rho2p=20*(-126-z1-5*z2-15*z3-35*z4+126*gammap)*t^3/tf^5+...
    30*(420+5*z1+24*z2+70*z3+160*z4-420*gammap)*t^4/tf^6+...
    42*(-540-10*z1-45*z2-126*z3-280*z4+540*gammap)*t^5/tf^7+...
    56*(315+10*z1+40*z2+105*z3+224*z4-315*gammap)*t^6/tf^8+...
    72*(-70-5*z1-15*z2-35*z3-70*z4+70*gammap)*t^7/tf^9+...
    90*z1*t^8/tf^10+110*z2*t^9/tf^11+132*z3*t^10/tf^12+156*z4*t^11/tf^13;
rho3p=60*(-126-z1-5*z2-15*z3-35*z4+126*gammap)*t^2/tf^5+...
    120*(420+5*z1+24*z2+70*z3+160*z4-420*gammap)*t^3/tf^6+...
    210*(-540-10*z1-45*z2-126*z3-280*z4+540*gammap)*t^4/tf^7+...
    336*(315+10*z1+40*z2+105*z3+224*z4-315*gammap)*t^5/tf^8+...
    504*(-70-5*z1-15*z2-35*z3-70*z4+70*gammap)*t^6/tf^9+...
    720*z1*t^7/tf^10+990*z2*t^8/tf^11+1320*z3*t^9/tf^12+1716*z4*t^10/tf^13;
rho4p=120*(-126-z1-5*z2-15*z3-35*z4+126*gammap)*t^1/tf^5+...
    360*(420+5*z1+24*z2+70*z3+160*z4-420*gammap)*t^2/tf^6+...
    840*(-540-10*z1-45*z2-126*z3-280*z4+540*gammap)*t^3/tf^7+...
    1680*(315+10*z1+40*z2+105*z3+224*z4-315*gammap)*t^4/tf^8+...
    3024*(-70-5*z1-15*z2-35*z3-70*z4+70*gammap)*t^5/tf^9+...
    5040*z1*t^6/tf^10+7920*z2*t^7/tf^11+11880*z3*t^8/tf^12+17160*z4*t^9/tf^13;
rhom=126*(gammam-1)*t^5/tf^5-420*(gammam-1)*t^6/tf^6+...
    540*(gammam-1)*t^7/tf^7-315*(gammam-1)*t^8/tf^8+70*(gammam-1)*t^9/tf^9+1;
rho1m=630*(gammam-1)*t^4/tf^5-2520*(gammam-1)*t^5/tf^6+...
    3780*(gammam-1)*t^6/tf^7-2520*(gammam-1)*t^7/tf^8+630*(gammam-1)*t^8/tf^9;
rho2m=2520*(gammam-1)*t^3/tf^5-12600*(gammam-1)*t^4/tf^6+...
    22680*(gammam-1)*t^5/tf^7-17640*(gammam-1)*t^6/tf^8+5040*(gammam-1)*t^7/tf^9;
rho3m=7560*(gammam-1)*t^2/tf^5-50400*(gammam-1)*t^3/tf^6+...
    113400*(gammam-1)*t^4/tf^7-105840*(gammam-1)*t^5/tf^8+35280*(gammam-1)*t^6/tf^9;
rho4m=15120*(gammam-1)*t^1/tf^5-151200*(gammam-1)*t^2/tf^6+...
    453600*(gammam-1)*t^3/tf^7-529200*(gammam-1)*t^4/tf^8+211680*(gammam-1)*t^5/tf^9;
wp=sqrt((sqrt(3)*w0)^2/rhop^4-rho2p/rhop);
w1p=1/2/wp*(-4*(sqrt(3)*w0)^2*rho1p/rhop^5-(rho3p*rhop-rho2p*rho1p)/rhop^2);
w2p=-1/4/wp^3*(-4*(sqrt(3)*w0)^2*rho1p/rhop^5-(rho3p*rhop-rho2p*rho1p)/rhop^2)^2+...
    1/(2*wp)*(-(4*(sqrt(3)*w0)^2*rho2p*rhop-20*(sqrt(3)*w0)^2*rho1p^2)/rhop^6-(rho4p*rhop^2-rho2p^2*rhop-2*rho3p*rho1p*rhop-2*rho2p*rho1p^2)/rhop^3);
wm=sqrt(w0^2/rhom^4-rho2m/rhom);
w1m=1/2/wm*(-4*w0^2*rho1m/rhom^5-(rho3m*rhom-rho2m*rho1m)/rhom^2);
w2m=-1/4/wm^3*(-4*w0^2*rho1m/rhom^5-(rho3m*rhom-rho2m*rho1m)/rhom^2)^2+...
    1/(2*wm)*(-(4*w0^2*rho2m*rhom-20*w0^2*rho1m^2)/rhom^6-(rho4m*rhom^2-rho2m^2*rhom-2*rho3m*rho1m*rhom-2*rho2m*rho1m^2)/rhom^3);
ddd=(4*2^(2/3)*Cc^(1/3)*(-2*m*wm*w1m+2*m*wp*w1p)^2)/(9*(-m*wm^2+m*wp^2)^(7/3))-...
    (2^(2/3)*Cc^(1/3)*(-2*m*w1m^2+2*m*w1p^2-2*m*wm*w2m+2*m*wp*w2p))/(3*(-m*wm^2+m*wp^2)^(4/3));
dydx=[y(2); -sqrt(m/2)*ddd-wp^2*y(1)];
end

[Jesús Parejo]

the code works fine, thank you! I have tried to write it with ode45 but MatLab takes a lot of time "busy". I don't really know the accuracy of the method because the solution oscillates roughly.

[Torsten]

From the description of your problem,

dydx=[y(2); -sqrt(m/2)*ddd-wp^2*y(1)];

should be

dydx=[y(2); -sqrt(m/2)*ddd+wp^2*y(1)];

shouldn't it ?

[Jesús Parejo]

No, no. Calculations are correct. I have mislead the interpretation.

Answer 2

[John D'Errico]

You have a simple classical ODE, with two INITIAL conditions, not boundary conditions at different ends. So use a tool like ODE45, NOT a boundary value solver. This is exctly what the ODE solvers (ODE45, etc.) are designed to solve.