Error solving bvp4c - Singular jacobian

2022 5 18
2 个回答
回答已采纳
更新时间：2022 5 18
12 次查看（30 天）

0 个投票

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

0 个投票

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更早的评论 隐藏 1更早的评论

Jesús Parejo
Jesús Parejo 2022-5-18
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
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

0 个投票

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.

类别

在 帮助中心 和 File Exchange 中查找有关 Boundary Value Problems 的更多信息

产品

版本

R2017b

标签

Community Treasure Hunt

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

Start Hunting!

Translated by