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 件のコメント
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 ExchangeOrdinary Differential Equations についてさらに検索

製品


リリース

R2017b

Community Treasure Hunt

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

Start Hunting!

Translated by