# Error solving bvp4c - Singular jacobian

3 views (last 30 days)
Jesús Parejo on 18 May 2022
Commented: Jesús Parejo on 18 May 2022
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?
%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
##### 2 CommentsShowHide 1 older comment
Jesús Parejo on 18 May 2022
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));

Torsten on 18 May 2022
Edited: Torsten on 18 May 2022
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 on 18 May 2022
No, no. Calculations are correct. I have mislead the interpretation.

John D'Errico on 18 May 2022
Edited: John D'Errico on 18 May 2022
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.

R2017b

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by