Error solving bvp4c - Singular jacobian

1 次查看（过去 30 天）
Jesús Parejo
Jesús Parejo 2022-5-18
评论： Jesús Parejo 2022-5-18
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 个评论
Torsten
Torsten 2022-5-18
Neither "solinit" nor "options" is supplied.
Jesús Parejo
Jesús Parejo 2022-5-18
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
Torsten 2022-5-18
编辑：Torsten 2022-5-18
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 个评论
显示 1更早的评论
Torsten
Torsten 2022-5-18
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 2022-5-18
No, no. Calculations are correct. I have mislead the interpretation.

请先登录，再进行评论。

更多回答（1 个）

John D'Errico
John D'Errico 2022-5-18
编辑：John D'Errico 2022-5-18
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.

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

产品


版本

R2017b

Community Treasure Hunt

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

Start Hunting!

Translated by