Error solving bvp4c - Singular jacobian

1 view (last 30 days)
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
  2 Comments
Torsten
Torsten on 18 May 2022
Neither "solinit" nor "options" is supplied.
Jesús Parejo
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));

Sign in to comment.

Accepted Answer

Torsten
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
  3 Comments
Torsten
Torsten on 18 May 2022
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
Jesús Parejo on 18 May 2022
No, no. Calculations are correct. I have mislead the interpretation.

Sign in to comment.

More Answers (1)

John D'Errico
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.

Products


Release

R2017b

Community Treasure Hunt

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

Start Hunting!