What is the desired solution?
Getting complex solution with Matlab ode45
13 次查看（过去 30 天）
显示 更早的评论
I'm trying to solve such an ODE:
with boundary condition
I have such a problem that if the solution of y would come with complex part such as 1.1462 + 0.0000i. However, if I add an abs in sqrt there will be no longer complex parts, but the solution is not what I desires. Is there anything I can do to avoid this?
My matlab code is
clearvars;
A=-4.32;
B=4;
C=1.2;
D=0.8;
a=0.2778;
odefun_mon_pc1=@(t,y) (t*D*y+A*t+C+sqrt(B*y))/((t*t-t)*D);
tspan2=[a 0.99];
options=odeset('Reltol',1e-10);
[t,y]=ode45(odefun_mon_pc1,tspan2,0.0001,options);
h_pc_mon=plot(t,y,'--blue','linewidth',1);
10 个评论
Kaiwei
2024-1-6
The desired solution should be a monotonically increasing one starts from a, and smoothly intersects with horizontal line at a.
Kaiwei
2024-1-6
Actually at around a, the numerator is also negative since A*a+C is negative. So around a, y' should be positive.
Paul
2024-1-6
Let's check
A=-4.32;
B=4;
C=1.2;
D=0.8;
a=0.2778;
odefun_mon_pc1=@(t,y) (t*D*y+A*t+C+sqrt(B*y))/((t*t-t)*D);
%tspan2=[a 0.99];
%options=odeset('Reltol',1e-10);
%[t,y]=ode45(odefun_mon_pc1,tspan2,0.0001,options);
%h_pc_mon=plot(t,y,'--blue','linewidth',1);
The derivative at t = a and y = 0 is
odefun_mon_pc1(a,0)
ans = 5.9812e-04
But, y just a bit more than 0, like the initial condition in the code
odefun_mon_pc1(a,0.0001)
ans = -0.1241
So it doesn't take much to make the derivative negative.
Kaiwei
2024-1-6
True, but when y is 0.0001 t should be larger than a, and since the y' at a is quite small, we expect a larger t so I guess y' should be positive anyway. For me it seems like a numerical error?
Kaiwei
2024-1-6
I actually want y(a) to be 0, since I feel it's possible to cause problems in sqrt, I use 0.0001 instead.
回答（2 个）
Paul
2024-1-6
编辑：Paul
2024-1-7
It looks like your differential equation is very sensitive (I mean that loosely)
A=-4.32;
B=4;
C=1.2;
D=0.8;
a=0.2778;
odefun_mon_pc1=@(t,y) (t*D*y+A*t+C+sqrt(B*y))/((t*t-t)*D);
tspan2=[a 0.99];
Here's the initial code:
options=odeset('Reltol',1e-10);
[t,y] = ode45(odefun_mon_pc1,tspan2,0.0001,options);
figure
h_pc_mon=plot(t,y,'--blue','linewidth',1);
Warning: Imaginary parts of complex X and/or Y arguments ignored.
Use the stated boundary condition:
options=odeset('Reltol',1e-10);
[t,y] = ode45(odefun_mon_pc1,tspan2,0*0.0001,options);
figure
h_pc_mon=plot(t,y,'--blue','linewidth',1);
Warning: Imaginary parts of complex X and/or Y arguments ignored.
Use the stated boundary condition with a smaller step size
options=odeset('Reltol',1e-10,'MaxStep',1e-6);
[t,y] = ode45(odefun_mon_pc1,tspan2,0*0.0001,options);
figure
h_pc_mon=plot(t,y,'--blue','linewidth',1);
If you cut tspan down to focus on the what's happening for t < 0.3, you'll probably find some interesting behavior.
I didn't try setting the MaxInitialStep to a small value; that's usually a good idea. And then maybe the MaxStep won't have to be set so small. You'll have to do some more analysis.
15 个评论
Kaiwei
2024-1-6
actually can you please check whether you get non-complex solutions? I find even though there is no warning that the complex part is ignored, but actually you can check y still brings 0.000i parts. What does that mean?
Walter Roberson
2024-1-7
编辑：Walter Roberson
2024-1-7
format long g
A=-4.32;
B=4;
C=1.2;
D=0.8;
a=0.2778;
odefun_mon_pc1=@(t,y) (t*D*y+A*t+C+sqrt(B*y))/((t*t-t)*D);
tspan2=[a 0.99];
options=odeset('Reltol',1e-10);
[t,y]=ode45(odefun_mon_pc1,tspan2,0.0001,options);
subplot(2,1,1)
plot(t,real(y),'--blue'); title('real');
subplot(2,1,2)
plot(t,imag(y),'--red'); title('imag')
[min(imag(y)), max(imag(y))]
ans = 1×2
1.0e+00 *
-1.32226601983809e-08 1.44587613035755e-06
With the default "format short", values in that range show up as 0.0000i
Torsten
2024-1-7
编辑：Torsten
2024-1-7
Maybe of interest: For given value of a, y0 must be chosen smaller or equal than the curve value in order to get an increasing solution for the differential equation.
syms y t real
A=-4.32;
B=4;
C=1.2;
D=0.8;
a=0.2778;
sol = solve(t*D*y+A*t+C+sqrt(B*y)==0,'ReturnConditions',1)
sol = struct with fields:
y: [2×1 sym]
parameters: [1×0 sym]
conditions: [2×1 sym]
sol.y
ans =
sol.conditions
ans =
y = sol.y(1)
y =
tnum = a:0.01:0.99;
format long
ynum = arrayfun(@(tnum)double(subs(y,t,tnum)),tnum);
ynum= ynum(:)
ynum = 72×1
0.000000002303975
0.000466314564157
0.001851362801389
0.004139423813455
0.007313383515290
0.011354852810903
0.016244285595462
0.021961097699154
0.028483785659285
0.035790044293692
plot(tnum,ynum)
xlabel('tstart')
ylabel('Frontier for y0')
Torsten
2024-1-7
编辑：Torsten
2024-1-7
For solutions with absolute value <1, it is important to restrict AbsTol, not RelTol in the computations:
A=-4.32;
B=4;
C=1.2;
D=0.8;
a=0.2778;
odefun_mon_pc1=@(t,y) (t*D*y+A*t+C+sqrt(B*y))/((t*t-t)*D);
tspan2=[a 0.99];
options=odeset('AbsTol',1e-10);
[t,y]=ode45(odefun_mon_pc1,tspan2,0,options);
subplot(2,1,1)
plot(t,real(y),'--blue'); title('real');
subplot(2,1,2)
plot(t,imag(y),'--red'); title('imag')
[min(imag(y)), max(imag(y))]
ans = 1×2
0 0
Paul
2024-1-7
In response to this comment, I don't see the solution have any element with a non-zero imaginary part.
A=-4.32;
B=4;
C=1.2;
D=0.8;
a=0.2778;
odefun_mon_pc1=@(t,y) (t*D*y+A*t+C+sqrt(B*y))/((t*t-t)*D);
tspan2=[a 0.99];
options=odeset('Reltol',1e-10,'MaxStep',1e-6);
[t,y] = ode45(odefun_mon_pc1,tspan2,0,options);
any(imag(y))
ans = logical
0
In response to this comment
sol = ode45(odefun_mon_pc1,tspan2,0,options);
f = @(t) deval(sol,t);
tval = 0.3:.01:.4;
[yval,yprime] = f(tval);
figure
plot(tval,yval,tval,yprime)
legend('y','yprime','Location','Northwest')
Kaiwei
2024-1-7
Thanks. Actually even setting Abstol you get the same result, the y still contains xx+0.000i part. I'm a little concerned that the appearance of 0.000i is a bug, there must be some cases there is negative numbers in sqrt. Do you know how to solve an ODE which can not be explicitly expressed as y'=@(t,y) ? (i.e. I want to squre on both side to get rid of the sqrt part. But after doing so I will get an ODE expressed as f(t,y,y')=0 ) I don't know how to solve that. Maybe you know how to do that?
Torsten
2024-1-7
编辑：Torsten
2024-1-7
Actually even setting Abstol you get the same result, the y still contains xx+0.000i part. I'm a little concerned that the appearance of 0.000i is a bug, there must be some cases there is negative numbers in sqrt.
I think once your function returned complex values (maybe in some intermediate step of the integration), the complex part in the solution remains - even if the imaginary part of the integration step didn't contribute to the final solution.
Do you know how to solve an ODE which can not be explicitly expressed as y'=@(t,y) ? (i.e. I want to squre on both side to get rid of the sqrt part. But after doing so I will get an ODE expressed as f(t,y,y')=0 ) I don't know how to solve that. Maybe you know how to do that?
Solve
z1' = z2
0 = f(t,z1,z2)
using ode15s with z1 = y and z2 = y'. Use the mass matrix option to do so.
But note that this squaring will complicate things for the solver because you will have to tell him which path to take when solving x^2 = c and c becomes 0 - the path of the positive or the negative square root.
Paul
2024-1-7
Please provide a specific example that illustrates the "appearance of 0.000i". As shown, when the code runs here on Answers as in the first part of this comment the imaginary part of all elements of y is identically zero.
Kaiwei
2024-1-7
I know the imaginary part is 0. But why do they appear? I'm afraid the reason is that actually it's very small but non zero. like 1e-9, so matlab treats it as 0.000 i. I guess if there is no bug, we should directly get real numbers.
Paul
2024-1-7
Without seeing an actual example what you're seeing that "appears," I'm afraid I can't be of any further help. Maybe someone else can. Why not post the actual code and output so we all know what you're talking about?
Kaiwei
2024-1-7
Sorry, I'm just referrring to the very original codes. Buy running this and you can see y contains the almost 0 imaginary part.
A=-4.32;
B=4;
C=1.2;
D=0.8;
a=5/18;
odefun_mon_pc1=@(t,y) (t*D*y+A*t+C+sqrt(B*y))/((t*t-t)*D);
tspan2=[a 0.99];
options=odeset('AbsTol',1e-10);
[t,y]=ode45(odefun_mon_pc1,tspan2,0,options);
y
Torsten
2024-1-7
Here is the "protocol" of the computation:
format long
main()
y =
0
dy =
1.383508692225195e-15
y =
1.970731270480778e-17
dy =
0.372052035790209
y =
0.005962133873538
dy =
-0.378907739937008
y =
-0.194879796143495
dy =
1.674569751746179 - 4.955806801539233i
y =
-0.607040471461461 + 0.102645099508424i
dy =
1.712248425195752 - 8.853224859229034i
y =
-0.525565056703810 + 0.074205897064859i
dy =
1.938343899800220 - 8.110876389255505i
y =
0.044288699548717 - 0.102168740039767i
dy =
-1.444861094730993 + 2.171921662465315i
y =
6.263417519238702e-18
dy =
0.120660369620546
y =
6.145350310098734e-04
dy =
-0.125149591382752
y =
-0.020269208618011
dy =
0.498156865621621 - 1.708403783043071i
y =
-0.062777574498281 + 0.011246006548785i
dy =
0.341710525171683 - 3.022754778844634i
y =
-0.053589014421031 + 0.007949372524639i
dy =
0.454526384023613 - 2.772586909505056i
y =
0.004922447429041 - 0.011337343018380i
dy =
-0.531195415656363 + 0.741600692506028i
y =
3.131708759619351e-18
dy =
0.060626617111502
y =
1.543886370637157e-04
dy =
-0.063189650213710
y =
-0.005104570386446
dy =
0.246177258146285 - 0.873172402330348i
y =
-0.015792105596448 + 0.002873940766314i
dy =
0.148309912877591 - 1.542930198052416i
y =
-0.013434669782905 + 0.002025248432693i
dy =
0.210303227142975 - 1.416765947698542i
y =
0.001263231849470 - 0.002904158325608i
dy =
-0.276910941487013 + 0.379385731866525i
y =
1.565854379809676e-18
dy =
0.030388541880217
y =
3.869295193076423e-05
dy =
-0.031753422780403
y =
-0.001280923169895
dy =
0.122445174153461 - 0.441614252687739i
y =
-0.003960720409051 + 0.000726759801614i
dy =
0.068571933330364 - 0.779955334290608i
y =
-0.003363612720182 + 0.000511532475800i
dy =
0.101026257547528 - 0.716654601507805i
y =
3.201622476449140e-04 - 7.352049543357601e-04i
dy =
-0.141541735028551 + 0.192062301447560i
y =
7.829271899048378e-19
dy =
0.015213220059255
y =
9.685301696691987e-06
dy =
-0.015916981842596
y =
-3.208360992182738e-04
dy =
0.061072270357478 - 0.222101698910611i
y =
-9.917974122841348e-04 + 1.827551824425714e-04i
dy =
0.032896057998398 - 0.392179726685343i
y =
-8.415374921489713e-04 + 1.285676070526490e-04i
dy =
0.049495488477978 - 0.360481754815441i
y =
8.060322930127133e-05 - 1.849756503688479e-04i
dy =
-0.071577114276770 + 0.096653242372683i
y =
3.914635949524189e-19
dy =
0.007611363809803
y =
2.422838640803803e-06
dy =
-0.007968635430208
y =
-8.028516210107499e-05
dy =
0.030499895984965 - 0.111379303593596i
y =
-2.481534074308239e-04 + 4.582392896680751e-05i
dy =
0.016101230182035 - 0.196650585781433i
y =
-2.104646168027045e-04 + 3.222951792616075e-05i
dy =
0.024495006224752 - 0.180790623032159i
y =
2.022232647710552e-05 - 4.639255219712682e-05i
dy =
-0.035994617795240 + 0.048485982471690i
y =
1.957317974762094e-19
dy =
0.003806871560430
y =
6.058990049396406e-07
dy =
-0.003986862411036
y =
-2.008083243529174e-05
dy =
0.015241057104542 - 0.055772367836485i
y =
-6.206396286788336e-05 + 1.147299785324150e-05i
dy =
0.007963996914436 - 0.098466756436168i
y =
-5.262623037663364e-05 + 8.068461078939240e-06i
dy =
0.012184524947846 - 0.090534149403213i
y =
5.064597100757451e-06 - 1.161683273272731e-05i
dy =
-0.018049396102513 + 0.024283316822590i
y =
9.786589873810472e-20
dy =
0.001903732738183
y =
1.514983830457891e-07
dy =
-0.001994067679724
y =
-5.021400623118931e-06
dy =
0.007618327035272 - 0.027906933817957i
y =
-1.551919344750512e-05 + 2.870383706163192e-06i
dy =
0.003960353947196 - 0.049268897484744i
y =
-1.315781520990902e-05 + 2.018508327617018e-06i
dy =
0.006076543469314 - 0.045301952604281i
y =
1.267281084102568e-06 - 2.906551654664398e-06i
dy =
-0.009037791259566 + 0.012151799069085i
y =
4.893294936905236e-20
dy =
9.519401219394521e-04
y =
3.787753037433204e-08
dy =
-9.971924487971494e-04
y =
-1.255498578340883e-06
dy =
0.003808615795127 - 0.013958659786063i
y =
-3.880196956626868e-06 + 7.178629847183327e-07i
dy =
0.001974767124185 - 0.024643342032173i
y =
-3.289609207000900e-06 + 5.048015998912719e-07i
dy =
0.003034345096395 - 0.022659711175699i
y =
3.169620927279541e-07 - 7.269310580455586e-07i
dy =
-0.004522177345437 + 0.006078440782650i
y =
2.446647468452618e-20
dy =
4.759881345279746e-04
y =
9.469742165442881e-09
dy =
-4.986354259235615e-04
y =
-3.138929353323804e-07
dy =
0.001904171299541 - 0.006980626320959i
y =
-9.700983066989318e-07 + 1.794990823893516e-07i
dy =
0.000986031284435 - 0.012323891944995i
y =
-8.224210217913041e-07 + 1.262222009289804e-07i
dy =
0.001516191227812 - 0.011332037176365i
y =
7.925828660738075e-08 - 1.817693480039817e-07i
dy =
-0.002261909763136 + 0.003039855506458i
y =
1.223323734226309e-20
dy =
2.379983252358120e-04
y =
2.367477897350908e-09
dy =
-2.493271833416031e-04
y =
-7.847542531909032e-08
dy =
0.000952051541917 - 0.003490635298130i
y =
-2.425304338343176e-07 + 4.487891230396562e-08i
dy =
0.000492677613895 - 0.006162498019250i
y =
-2.056073416401360e-07 + 3.155826106954752e-08i
dy =
0.000757850291163 - 0.005666560848646i
y =
1.981679929002379e-08 - 4.544688049533326e-08i
dy =
-0.001131159959216 + 0.001520085829139i
y =
6.116618671131545e-21
dy =
1.190000426484157e-04
y =
5.918738513713144e-10
dy =
-1.246657245091535e-04
y =
-1.961909659116244e-08
dy =
0.000476017249224 - 0.001745396703070i
y =
-6.063324237884449e-08 + 1.122023627255852e-08i
dy =
0.000246254317711 - 0.003081384582084i
y =
-5.140200617536733e-08 + 7.889898484132458e-09i
dy =
0.000378863791010 - 0.002833413567497i
y =
4.954480487493706e-09 - 1.136227712690536e-08i
dy =
-5.656310282254856e-04 + 7.600818660153602e-04i
y =
3.058309335565772e-21
dy =
5.950011085744896e-05
y =
1.479686855003989e-10
dy =
-6.233322929042293e-05
y =
-4.904792266231969e-09
dy =
2.380064948384422e-04 - 8.727170572473672e-04i
y =
-1.515835507313348e-08 + 2.805119192727366e-09i
dy =
0.000123106047744 - 0.001540724436017i
y =
-1.285049104964707e-08 + 1.972514243038194e-09i
dy =
0.000189416536322 - 0.001416738331473i
y =
1.238655965218389e-09 - 2.840635008699860e-09i
dy =
-2.828281097695261e-04 + 3.800502607453933e-04i
y =
1.529154667782886e-21
dy =
2.974998554798552e-05
y =
3.699208448343836e-11
dy =
-3.116658887906363e-05
y =
-1.226196118504894e-09
dy =
1.190027150621082e-04 - 4.363624557087092e-04i
y =
-3.789582035860693e-09 + 7.012861094722451e-10i
dy =
6.154775382903659e-05 - 7.703690152113238e-04i
y =
-3.212610192462292e-09 + 4.931327934659746e-10i
dy =
9.470441337666793e-05 - 7.083758204797021e-04i
y =
3.096687391795505e-10 - 7.101655877030202e-10i
dy =
-1.414170847732855e-04 + 1.900271712024422e-04i
y =
1.533642674400130e-10 - 5.001337046489064e-10i
dy =
-5.068080282980875e-05 + 1.694374462745171e-04i
y =
1.880364294813022e-10 - 4.207196862819207e-10i
dy =
-3.107853169619963e-05 + 1.455349260455315e-04i
y =
-1.952352244550604e-11 - 3.194770716137021e-10i
dy =
1.149906276898515e-04 + 1.623790275170763e-04i
y =
-6.219946268655014e-10 - 8.280704800216605e-10i
dy =
1.017594622116897e-04 + 3.587518707200944e-04i
y =
-4.083600624752463e-10 - 9.234447260919244e-10i
dy =
8.142172857915304e-05 + 3.318150495992067e-04i
y =
4.526360649636436e-10 - 6.792365410744312e-11i
dy =
3.164068063595222e-05 + 1.983660212809084e-05i
y =
4.968760198415248e-10 - 4.018814509060903e-11i
dy =
5.713509847636958e-05 + 1.122409970957442e-05i
y =
5.590980385847979e-10 - 3.986762785910936e-11i
dy =
5.911254087223089e-05 + 1.049861967195837e-05i
y =
6.470633371323969e-10 + 3.568921347722917e-11i
dy =
1.309304882661696e-04 - 8.738295526763915e-06i
y =
2.672054766348865e-10 + 1.703655906023332e-10i
dy =
2.518016286169323e-04 - 6.211254990232522e-05i
y =
2.392767173857214e-10 + 2.381306722973853e-10i
dy =
2.740338440321297e-04 - 8.736202474152284e-05i
y =
9.377491449023052e-10 - 2.077368415140676e-12i
dy =
1.040699535173227e-04 + 4.226658719070449e-07i
y =
1.131931268053030e-09 - 1.288724300135851e-12i
dy =
1.166401436486308e-04 + 2.386578249497965e-07i
y =
1.255408614923137e-09 - 1.280656449835115e-12i
dy =
1.194757525773124e-04 + 2.251982787891306e-07i
y =
1.787689803370570e-09 + 9.359270766798996e-13i
dy =
1.596777182009342e-04 - 1.379167846705754e-07i
y =
1.701860482045778e-09 + 4.758676962501275e-12i
dy =
1.948005317660836e-04 - 7.186941903966038e-07i
y =
1.840087150086763e-09 + 5.381890287541470e-12i
dy =
2.022318289762797e-04 - 7.816901731783460e-07i
y =
2.158025114786880e-09 - 4.052817816603760e-13i
dy =
1.578926058343788e-04 + 5.435614502240691e-08i
y =
2.716234230729141e-09 - 2.131125822311487e-13i
dy =
1.824811205844441e-04 + 2.547665723461021e-08i
y =
3.093134531790896e-09 - 2.318901853373331e-13i
dy =
1.864641755727737e-04 + 2.597757228921311e-08i
y =
4.563930095362158e-09 + 4.856333788727726e-13i
dy =
2.755515348996356e-04 - 4.478651681369425e-08i
y =
3.955139642687552e-09 + 1.950497105815681e-12i
dy =
3.759932983358271e-04 - 1.932282741881028e-07i
y =
4.296045329305949e-09 + 2.288631616413395e-12i
dy =
3.957796700454491e-04 - 2.175434973018268e-07i
y =
5.837821898554267e-09 - 2.927958145685065e-14i
dy =
2.604411132178439e-04 + 2.387509717550282e-09i
y =
7.015789542074857e-09 - 1.848094386420300e-14i
dy =
2.905275378158513e-04 + 1.374636703140663e-09i
y =
7.757863398387546e-09 - 1.823546992797686e-14i
dy =
2.975822711275659e-04 + 1.289870770004645e-09i
y =
1.099599209508535e-08 + 1.117125159425661e-14i
dy =
3.927756622486992e-04 - 6.637058943682332e-10i
y =
1.056380081721246e-08 + 6.056085475219572e-14i
dy =
4.728093571104030e-04 - 3.670885270199564e-09i
y =
1.140860157977535e-08 + 6.830817900799865e-14i
dy =
4.902099993699877e-04 - 3.984220730044149e-09i
y =
1.318556707868458e-08 - 6.062471291534953e-15i
dy =
3.903137586007077e-04 + 3.289182360150243e-10i
y =
1.558039589650263e-08 - 4.044343944093738e-15i
dy =
4.309227854862468e-04 + 2.018559549034899e-10i
y =
1.705811845217644e-08 - 3.912341241138842e-15i
dy =
4.413990048913700e-04 + 1.866171279406365e-10i
y =
2.368613179414105e-08 + 1.040800248216085e-15i
dy =
5.639246147271665e-04 - 4.212967926019293e-11i
y =
2.322867364824396e-08 + 8.535940083533973e-15i
dy =
6.559142672277084e-04 - 3.489026130972983e-10i
y =
2.497002453165162e-08 + 9.599041640615555e-15i
dy =
6.777506602729647e-04 - 3.784257895407139e-10i
y =
2.785922920630100e-08 - 1.481942161835148e-15i
dy =
5.669655527584294e-04 + 5.531079475239476e-11i
y =
3.228084084586639e-08 - 1.050588170960497e-15i
dy =
6.178729961737999e-04 + 3.642646390086110e-11i
y =
3.493828678541850e-08 - 1.000594211783202e-15i
dy =
6.324836041593436e-04 + 3.334732956591502e-11i
y =
4.721844953293655e-08 - 5.251995556896841e-17i
dy =
7.786078530353374e-04 + 1.505589432703325e-12i
y =
4.718785405821270e-08 + 1.171429287270277e-15i
dy =
8.727484352702047e-04 - 3.359198708562287e-11i
y =
5.040295304993534e-08 + 1.332711163108839e-15i
dy =
8.986491215644930e-04 - 3.697765051349517e-11i
y =
5.433764453625484e-08 - 4.295297449077720e-16i
dy =
7.915276587965475e-04 + 1.147823042814519e-11i
y =
6.203073741874014e-08 - 3.179694073390918e-16i
dy =
8.542650061542439e-04 + 7.952545519226411e-12i
y =
6.656326717900159e-08 - 3.007398170777042e-16i
dy =
8.736576902279902e-04 + 7.260961045302960e-12i
y =
8.791881820770416e-08 - 7.232583959546645e-17i
dy =
0.001047622226649 + 0.000000000001519i
y =
8.876196973263427e-08 + 1.806642602244033e-16i
dy =
0.001146168961779 - 0.000000000003777i
y =
9.430360174700384e-08 + 2.141514210816973e-16i
dy =
0.001177334867269 - 0.000000000004344i
y =
9.959779652033275e-08 - 1.405718604017924e-16i
dy =
0.001071419768284 + 0.000000000002774i
y =
1.122733652517125e-07 - 1.077491311612818e-16i
dy =
0.001147021989341 + 0.000000000002003i
y =
1.196173740915352e-07 - 1.016058253473297e-16i
dy =
0.001171806627979 + 0.000000000001830i
y =
1.547177914918584e-07 - 3.757729233055981e-17i
dy =
0.001375957686682 + 0.000000000000595i
y =
1.571675428903223e-07 + 2.311747542081756e-17i
dy =
0.001478761577743 - 0.000000000000363i
y =
1.661791717131384e-07 + 3.169210026289701e-17i
dy =
0.001515962106241 - 0.000000000000484i
y =
1.730471056647400e-07 - 5.090057355869628e-17i
dy =
0.001412098247470 + 0.000000000000762i
y =
1.953279764305383e-07 - 3.887643455633899e-17i
dy =
0.001513044533853 + 0.000000000000548i
y =
2.082602967342317e-07 - 3.666723903467921e-17i
dy =
0.001545936870256 + 0.000000000000500i
y =
2.699813305988067e-07 - 1.310574319555488e-17i
dy =
0.001819214849423 + 0.000000000000157i
y =
2.741036512867933e-07 + 9.623451643950829e-18i
dy =
0.001958598027056 - 0.000000000000114i
y =
2.899697089998235e-07 + 1.280263918685518e-17i
dy =
0.002008286284741 - 0.000000000000148i
y =
3.023655456484799e-07 - 1.823743692445759e-17i
dy =
0.001866441392225 + 0.000000000000207i
y =
3.412454194595944e-07 - 1.393522733327800e-17i
dy =
0.001999670340534 + 0.000000000000149i
y =
3.638075627788081e-07 - 1.314305724203303e-17i
dy =
0.002043112235569 + 0.000000000000136i
y =
4.715059148129636e-07 - 4.717428047450815e-18i
dy =
0.002403652942016 + 0.000000000000043i
y =
4.787372215848366e-07 + 3.393865088988920e-18i
dy =
0.002587260368816 - 0.000000000000031i
y =
5.064182930011868e-07 + 4.529626307766129e-18i
dy =
0.002652826793885 - 0.000000000000040i
y =
5.279820937383298e-07 - 6.542876330080363e-18i
dy =
0.002466137501859 + 0.000000000000056i
y =
5.958886160960141e-07 - 4.999208393492945e-18i
dy =
0.002642200756944 + 0.000000000000040i
y =
6.352958812704481e-07 - 4.715053994476078e-18i
dy =
0.002699604997647 + 0.000000000000037i
y =
8.233994505256927e-07 - 1.691714331689609e-18i
dy =
0.003176054001929 + 0.000000000000012i
y =
8.360232915735843e-07 + 1.219521057971129e-18i
dy =
0.003418713120302 - 0.000000000000008i
y =
8.843712101250965e-07 + 1.627068266172843e-18i
dy =
0.003505349381809 - 0.000000000000011i
y =
9.220443580594382e-07 - 2.347140589559471e-18i
dy =
0.003258602372202 + 0.000000000000015i
y =
1.040644364069050e-06 - 1.793398284713334e-18i
dy =
0.003491232031364 + 0.000000000000011i
y =
1.109469497881838e-06 - 1.691471214669743e-18i
dy =
0.003567079947758 + 0.000000000000010i
y =
1.437994953083668e-06 - 6.069599909703393e-19i
dy =
0.004196561063946 + 0.000000000000003i
y =
1.460046915726589e-06 + 4.373312977636748e-19i
dy =
0.004517127973924 - 0.000000000000002i
y =
1.544487017101725e-06 + 5.835057252878616e-19i
dy =
0.004631580554002 - 0.000000000000003i
y =
1.610275472384083e-06 - 8.420700904924294e-19i
dy =
0.004305630936470 + 0.000000000000004i
y =
1.817431951411863e-06 - 6.434118196959995e-19i
dy =
0.004612999065390 + 0.000000000000003i
y =
1.937647116576954e-06 - 6.068485461531244e-19i
dy =
0.004713214812273 + 0.000000000000003i
y =
2.511474781323582e-06 - 2.177776356562922e-19i
dy =
0.005544863156345 + 0.000000000000001i
y =
2.549998428965226e-06 + 1.568762614890002e-19i
dy =
0.005968340240279 - 0.000000000000001i
y =
2.697487892663802e-06 + 2.093090763254276e-19i
dy =
0.006119531775146 - 0.000000000000001i
y =
2.812386649625631e-06 - 3.021327816731008e-19i
dy =
0.005688967798223 + 0.000000000000001i
y =
3.174261656700916e-06 - 2.308570419999297e-19i
dy =
0.006095072993890 + 0.000000000000001i
y =
3.384260536834194e-06 - 2.177403386224525e-19i
dy =
0.006227481041714 + 0.000000000000001i
y =
4.386657410815355e-06 - 7.814919979577551e-20i
dy =
0.007326160030987 + 0.000000000000000i
y =
4.453967413350881e-06 + 5.627455492680053e-20i
dy =
0.007885530390488 - 0.000000000000000i
y =
4.711610672109094e-06 + 7.508281576759433e-20i
dy =
0.008085234113719 - 0.000000000000000i
y =
4.912292977196762e-06 - 1.084180932055545e-19i
dy =
0.007516544984755 + 0.000000000000000i
y =
5.544530479019582e-06 - 8.284244115666805e-20i
dy =
0.008053083490775 + 0.000000000000000i
y =
5.911420193989197e-06 - 7.813660966718626e-20i
dy =
0.008228017754643 + 0.000000000000000i
y =
7.662713564511907e-06 - 2.804845374443637e-20i
dy =
0.009679357709568 + 0.000000000000000i
y =
7.780344660347781e-06 + 2.018801115135131e-20i
dy =
0.010418140949827 - 0.000000000000000i
y =
8.230474053515427e-06 + 2.693509270000566e-20i
dy =
0.010681888932555 - 0.000000000000000i
y =
8.581018899156376e-06 - 3.891140999419384e-20i
dy =
0.009930871202460 + 0.000000000000000i
y =
9.685816631900941e-06 - 2.973279742852842e-20i
dy =
0.010639697827190 + 0.000000000000000i
y =
1.032692864195810e-05 - 2.804433933174525e-20i
dy =
0.010870803960216 + 0.000000000000000i
y =
1.338718723286769e-05 - 1.006907958998351e-20i
dy =
0.012787812263328 + 0.000000000000000i
y =
1.359281602430164e-05 + 7.242827921072521e-21i
dy =
0.013763399006427 - 0.000000000000000i
y =
1.437938297672545e-05 + 9.663365312243331e-21i
dy =
0.014111671050880 - 0.000000000000000i
y =
1.499177727477740e-05 - 1.396837304529752e-20i
dy =
0.013120054935463 + 0.000000000000000i
y =
1.692281785572957e-05 - 1.067369271637113e-20i
dy =
0.014056426505277 + 0.000000000000000i
y =
1.804338269043093e-05 - 1.006779416641831e-20i
dy =
0.014361719008736 + 0.000000000000000i
y =
2.339225605974263e-05 - 3.615738043200563e-21i
dy =
0.016893472159108 + 0.000000000000000i
y =
2.375184148787283e-05 + 2.598766278751253e-21i
dy =
0.018181497699951 - 0.000000000000000i
y =
2.512663837209126e-05 + 3.467217948367936e-21i
dy =
0.018641279084035 - 0.000000000000000i
y =
2.619664941617839e-05 - 5.015768267447482e-21i
dy =
0.017332309147307 + 0.000000000000000i
y =
2.957293304506768e-05 - 3.832829456684031e-21i
dy =
0.018569153612629 + 0.000000000000000i
y =
3.153212529455772e-05 - 3.615367661846114e-21i
dy =
0.018972407021990 + 0.000000000000000i
y =
4.088411122686798e-05 - 1.298880477094218e-21i
dy =
0.022315455844949 + 0.000000000000000i
y =
4.151321881441775e-05 + 9.325791450248254e-22i
dy =
0.024015517467022 - 0.000000000000000i
y =
4.391690634912977e-05 + 1.244202699256956e-21i
dy =
0.024622333594720 - 0.000000000000000i
y =
4.578688348550062e-05 - 1.801729773489379e-21i
dy =
0.022895002929960 + 0.000000000000000i
y =
5.169253329032936e-05 - 1.376856988334038e-21i
dy =
0.024528530825700 + 0.000000000000000i
y =
5.511938845140728e-05 - 1.298790842829217e-21i
dy =
0.025061107855570 + 0.000000000000000i
y =
7.147717213282131e-05 - 4.668272522585244e-22i
dy =
0.029474408911594 + 0.000000000000000i
y =
7.257849533557888e-05 + 3.347217951049519e-22i
dy =
0.031717515755960 - 0.000000000000000i
y =
7.678281208731770e-05 + 4.465581244454974e-22i
dy =
0.032518080166005 - 0.000000000000000i
y =
8.005172008190541e-05 - 6.475168188291440e-22i
dy =
0.030239635274729 + 0.000000000000000i
y =
9.038723073271329e-05 - 4.948487331738420e-22i
dy =
0.032396701149882 + 0.000000000000000i
y =
9.638439998011620e-05 - 4.668157204994411e-22i
dy =
0.033099949982845 + 0.000000000000000i
y =
1.250115117491645e-04 - 1.678886259790811e-22i
dy =
0.038924352602890 + 0.000000000000000i
y =
1.269410416996646e-04 + 1.201683118035633e-22i
dy =
0.041882591851017 - 0.000000000000000i
y =
1.342987647317007e-04 + 1.603129006064665e-22i
dy =
0.042938236134581 - 0.000000000000000i
y =
1.400151617819276e-04 - 2.328540284838309e-22i
dy =
0.039934480394832 + 0.000000000000000i
y =
1.581160674742267e-04 - 1.779646755521442e-22i
dy =
0.042782213008746 + 0.000000000000000i
y =
1.686186453096881e-04 - 1.678943361430487e-22i
dy =
0.043710601382799 + 0.000000000000000i
y =
2.187524227348342e-04 - 6.042877624015380e-23i
dy =
0.051394173551753 + 0.000000000000000i
y =
2.221364683368635e-04 + 4.315607398676881e-23i
dy =
0.055293184662989 - 0.000000000000000i
y =
2.350215797577061e-04 + 5.757041652441292e-23i
dy =
0.056684273885666 - 0.000000000000000i
y =
2.450224806973521e-04 - 8.380424159287387e-23i
dy =
0.052727096699245 + 0.000000000000000i
y =
2.767516101623589e-04 - 6.405494065913138e-23i
dy =
0.056485446065454 + 0.000000000000000i
y =
2.951605080027951e-04 - 6.043552803015446e-23i
dy =
0.057710640321346 + 0.000000000000000i
y =
3.830357646250116e-04 - 2.177315546310597e-23i
dy =
0.067841394575844 + 0.000000000000000i
y =
3.889786139576488e-04 + 1.550575823898543e-23i
dy =
0.072976358998487 - 0.000000000000000i
y =
4.115631435579221e-04 + 2.068338916851002e-23i
dy =
0.074807930947693 - 0.000000000000000i
y =
4.290699324557458e-04 - 3.019218357789487e-23i
dy =
0.069599380546142 + 0.000000000000000i
y =
4.847511715669524e-04 - 2.307954930901419e-23i
dy =
0.074557366315225 + 0.000000000000000i
y =
5.170541173325157e-04 - 2.177782167902278e-23i
dy =
0.076173504540351 + 0.000000000000000i
y =
6.712548933613435e-04 - 7.855437525118896e-24i
dy =
0.089521352431685 + 0.000000000000000i
y =
6.817088354575994e-04 + 5.574645876676866e-24i
dy =
0.096277322724043 - 0.000000000000000i
y =
7.213374620663292e-04 + 7.435441004473227e-24i
dy =
0.098686218323827 - 0.000000000000000i
y =
7.520056212654031e-04 - 1.089144737016721e-23i
dy =
0.091838316474623 + 0.000000000000000i
y =
8.498582541463439e-04 - 8.326753280321716e-24i
dy =
0.098374873696607 + 0.000000000000000i
y =
9.066197724779052e-04 - 7.858177912458822e-24i
dy =
0.100505311044593 + 0.000000000000000i
y =
0.001177580142772 - 0.000000000000000i
dy =
0.118075253690363 + 0.000000000000000i
y =
0.001196007857991 + 0.000000000000000i
dy =
0.126952759704512 - 0.000000000000000i
y =
0.001265637085649 + 0.000000000000000i
dy =
0.130116528113751 - 0.000000000000000i
y =
0.001319408650183 - 0.000000000000000i
dy =
0.121125871580884 + 0.000000000000000i
y =
0.001491666146368 - 0.000000000000000i
dy =
0.129736197558780 + 0.000000000000000i
y =
0.001591570583516 - 0.000000000000000i
dy =
0.132541988692257 + 0.000000000000000i
y =
0.002068491337367 - 0.000000000000000i
dy =
0.155640482876983 + 0.000000000000000i
y =
0.002101057449425 + 0.000000000000000i
dy =
0.167287482951030 - 0.000000000000000i
y =
0.002223595345533 + 0.000000000000000i
dy =
0.171435343327616 - 0.000000000000000i
y =
0.002317972993879 - 0.000000000000000i
dy =
0.159651261080661 + 0.000000000000000i
y =
0.002621812289381 - 0.000000000000000i
dy =
0.170979288568378 + 0.000000000000000i
y =
0.002797985659598 - 0.000000000000000i
dy =
0.174669548963467 + 0.000000000000000i
y =
0.003639021779244 - 0.000000000000000i
dy =
0.204985366818227 + 0.000000000000000i
y =
0.003696745605829 + 0.000000000000000i
dy =
0.220236856954924 - 0.000000000000000i
y =
0.003912785907145 + 0.000000000000000i
dy =
0.225662991101011 - 0.000000000000000i
y =
0.004078623143704 - 0.000000000000000i
dy =
0.210247360474903 + 0.000000000000000i
y =
0.004615695605331 - 0.000000000000000i
dy =
0.225124014416973 + 0.000000000000000i
y =
0.004926984198301 - 0.000000000000000i
dy =
0.229967920849140 + 0.000000000000000i
y =
0.006413109472827 - 0.000000000000000i
dy =
0.269667081575811 + 0.000000000000000i
y =
0.006515763154695 + 0.000000000000000i
dy =
0.289594356573050 - 0.000000000000000i
y =
0.006897355033774 + 0.000000000000000i
dy =
0.296673861459230 - 0.000000000000000i
y =
0.007189077325349 - 0.000000000000000i
dy =
0.276548914532650 + 0.000000000000000i
y =
0.008140256575540 - 0.000000000000000i
dy =
0.296033296848327 + 0.000000000000000i
y =
0.008691238948847 - 0.000000000000000i
dy =
0.302372537902084 + 0.000000000000000i
y =
0.011321843829659 - 0.000000000000000i
dy =
0.354203066226909 + 0.000000000000000i
y =
0.011504990932320 + 0.000000000000000i
dy =
0.380176651139166 - 0.000000000000000i
y =
0.012179984715132 + 0.000000000000000i
dy =
0.389384837531925 - 0.000000000000000i
y =
0.012693503788658 - 0.000000000000000i
dy =
0.363160607085906 + 0.000000000000000i
y =
0.014379788962956 - 0.000000000000000i
dy =
0.388575524387946 + 0.000000000000000i
y =
0.015355693466165 - 0.000000000000000i
dy =
0.396833888968967 + 0.000000000000000i
y =
0.020015462693770 - 0.000000000000000i
dy =
0.464228902743291 + 0.000000000000000i
y =
0.020342990971860 + 0.000000000000000i
dy =
0.498004979493091 - 0.000000000000000i
y =
0.021537254388306 + 0.000000000000000i
dy =
0.509942096782366 - 0.000000000000000i
y =
0.022440804377068 - 0.000000000000000i
dy =
0.475803101020528 + 0.000000000000000i
y =
0.025426219081275 - 0.000000000000000i
dy =
0.508744723116854 + 0.000000000000000i
y =
0.027151454243986 - 0.000000000000000i
dy =
0.519427201039156 + 0.000000000000000i
y =
0.035390327589322 - 0.000000000000000i
dy =
0.606583696514160 + 0.000000000000000i
y =
0.035975929218288 + 0.000000000000000i
dy =
0.650440094194655 - 0.000000000000000i
y =
0.038083309311274 + 0.000000000000000i
dy =
0.665868204284890 - 0.000000000000000i
y =
0.039668471200948 - 0.000000000000000i
dy =
0.621368205592725 + 0.000000000000000i
y =
0.044920057808670 - 0.000000000000000i
dy =
0.663650245951446 + 0.000000000000000i
y =
0.047947873245718 - 0.000000000000000i
dy =
0.677316421798375 + 0.000000000000000i
y =
0.062410534698806 - 0.000000000000000i
dy =
0.789219363312337 + 0.000000000000000i
y =
0.063451166861891 + 0.000000000000000i
dy =
0.846213174169114 - 0.000000000000000i
y =
0.067138174873697 + 0.000000000000000i
dy =
0.866126765806761 - 0.000000000000000i
y =
0.069897797453880 - 0.000000000000000i
dy =
0.807759666054246 + 0.000000000000000i
y =
0.078989966966805 - 0.000000000000000i
dy =
0.861233811848343 + 0.000000000000000i
y =
0.084213196820837 - 0.000000000000000i
dy =
0.878420479327049 + 0.000000000000000i
y =
0.109170624124542 - 0.000000000000000i
dy =
1.020816673464679 + 0.000000000000000i
y =
0.110986149805521 + 0.000000000000000i
dy =
1.095371041804911 - 0.000000000000000i
y =
0.117313209649450 + 0.000000000000000i
dy =
1.121172052833214 - 0.000000000000000i
y =
0.122043346039916 - 0.000000000000000i
dy =
1.043359949332595 + 0.000000000000000i
y =
0.136905428873742 - 0.000000000000000i
dy =
1.107751035939996 + 0.000000000000000i
y =
0.145368337453539 - 0.000000000000000i
dy =
1.128527184093873 + 0.000000000000000i
y =
0.185937613825935 - 0.000000000000000i
dy =
1.301938231632023 + 0.000000000000000i
y =
0.189145859760665 + 0.000000000000000i
dy =
1.395074123938578 - 0.000000000000000i
y =
0.199318677933949 + 0.000000000000000i
dy =
1.427901740275793 - 0.000000000000000i
y =
0.206579693068181 - 0.000000000000000i
dy =
1.328785162137657 + 0.000000000000000i
y =
0.225507499488853 - 0.000000000000000i
dy =
1.390262314231399 + 0.000000000000000i
y =
0.235956574061491 - 0.000000000000000i
dy =
1.411882054709986 + 0.000000000000000i
y =
0.286987377386788 - 0.000000000000000i
dy =
1.574416832359001 + 0.000000000000000i
y =
0.292978367156633 - 0.000000000000000i
dy =
1.648800146450195 + 0.000000000000000i
y =
0.305463776512641 - 0.000000000000000i
dy =
1.681251656236709 + 0.000000000000000i
y =
0.311206766591890 - 0.000000000000000i
dy =
1.609012152413721 + 0.000000000000000i
y =
0.334126250807384 - 0.000000000000000i
dy =
1.670269479000037 + 0.000000000000000i
y =
0.346567641573676 - 0.000000000000000i
dy =
1.692490033772320 + 0.000000000000000i
y =
0.407736092408705 - 0.000000000000000i
dy =
1.855883799010483 + 0.000000000000000i
y =
0.415768640369989 - 0.000000000000000i
dy =
1.928155671547469 + 0.000000000000000i
y =
0.430500904022383 - 0.000000000000000i
dy =
1.962364248708054 + 0.000000000000000i
y =
0.435890183975403 - 0.000000000000000i
dy =
1.893608124432500 + 0.000000000000000i
y =
0.462863579703430 - 0.000000000000000i
dy =
1.957400652059764 + 0.000000000000000i
y =
0.477372552822671 - 0.000000000000000i
dy =
1.980398427611395 + 0.000000000000000i
y =
0.548799874983569 - 0.000000000000000i
dy =
2.156051882715046 + 0.000000000000000i
y =
0.558370541792580 - 0.000000000000000i
dy =
2.240384606364275 + 0.000000000000000i
y =
0.575383094678526 - 0.000000000000000i
dy =
2.279717664030780 + 0.000000000000000i
y =
0.581342256470338 - 0.000000000000000i
dy =
2.194294047260958 + 0.000000000000000i
y =
0.612598756121322 - 0.000000000000000i
dy =
2.264548340043944 + 0.000000000000000i
y =
0.629352830988661 - 0.000000000000000i
dy =
2.288547498367116 + 0.000000000000000i
y =
0.711556130788430 - 0.000000000000000i
dy =
2.497009384422289 + 0.000000000000000i
y =
0.721909172436788 - 0.000000000000000i
dy =
2.624561760756301 + 0.000000000000000i
y =
0.741205965780401 - 0.000000000000000i
dy =
2.680335724079038 + 0.000000000000000i
y =
0.749331025942823 - 0.000000000000000i
dy =
2.529212054066199 + 0.000000000000000i
y =
0.785358246535188 - 0.000000000000000i
dy =
2.613933955875900 + 0.000000000000000i
y =
0.804729525307872 - 0.000000000000000i
dy =
2.638131409299489 + 0.000000000000000i
y =
0.898494814264489 - 0.000000000000000i
dy =
2.952181156160045 + 0.000000000000000i
y =
0.906922344882561 + 0.000000000000000i
dy =
3.288820829732756 - 0.000000000000000i
y =
0.927235110451703 + 0.000000000000000i
dy =
3.468027514816661 - 0.000000000000000i
y =
0.943879756234952 - 0.000000000000000i
dy =
2.926187983113604 + 0.000000000000000i
y =
0.981083739149486 - 0.000000000000000i
dy =
3.040157641411843 + 0.000000000000000i
y =
1.001315885914664 - 0.000000000000000i
dy =
3.052896241084051 + 0.000000000000000i
y =
1.094286960361290 + 0.000000000000000i
dy =
4.176378672346862 - 0.000000000000000i
y =
1.081228260216912 + 0.000000000000000i
dy =
8.676566775219545 - 0.000000000000000i
y =
1.045836017143651 + 0.000000000000000i
dy =
25.653458231737751 - 0.000000000000000i
y =
1.256612236625872 - 0.000000000000000i
dy =
-20.253943508377660 + 0.000000000000000i
y =
0.956794372763178 - 0.000000000000000i
dy =
2.959694562493763 + 0.000000000000000i
y =
0.963418046006117 - 0.000000000000000i
dy =
2.972775930415163 + 0.000000000000000i
y =
0.996433157601662 - 0.000000000000000i
dy =
3.064880409889199 + 0.000000000000000i
y =
1.001912707466349 - 0.000000000000000i
dy =
3.103140002326774 + 0.000000000000000i
y =
1.009439108431380 - 0.000000000000000i
dy =
3.125040007225063 + 0.000000000000000i
y =
1.010223046653542 - 0.000000000000000i
dy =
3.089857307674118 + 0.000000000000000i
y =
1.023860011361834 - 0.000000000000000i
dy =
3.129277314758268 + 0.000000000000000i
y =
1.030874219712640 - 0.000000000000000i
dy =
3.142753813020164 + 0.000000000000000i
y =
1.065673644442202 - 0.000000000000000i
dy =
3.266424403486700 + 0.000000000000000i
y =
1.071077676673559 - 0.000000000000000i
dy =
3.354563130607045 + 0.000000000000000i
y =
1.078933237120034 - 0.000000000000000i
dy =
3.397961639318861 + 0.000000000000000i
y =
1.080476188169096 - 0.000000000000000i
dy =
3.280575897194256 + 0.000000000000000i
y =
1.093228523438183 - 0.000000000000000i
dy =
3.325774388342657 + 0.000000000000000i
y =
1.099802349844296 - 0.000000000000000i
dy =
3.336664983199894 + 0.000000000000000i
y =
1.132081962200902 + 0.000000000000000i
dy =
3.553649773565003 - 0.000000000000000i
y =
1.136131521329597 + 0.000000000000000i
dy =
3.936566663465324 - 0.000000000000000i
y =
1.142487042185821 + 0.000000000000000i
dy =
4.319009191872722 - 0.000000000000000i
y =
1.146715761062576 - 0.000000000000000i
dy =
3.397072800615574 + 0.000000000000000i
function main
A=-4.32;
B=4;
C=1.2;
D=0.8;
a=5/18;
%odefun_mon_pc1=@(t,y) (t*D*y+A*t+C+sqrt(B*y))/((t*t-t)*D);
tspan2=[a 0.99];
options=odeset('AbsTol',1e-10);
[t,y]=ode45(@odefun_mon_pc1,tspan2,0,options);
function dy = odefun_mon_pc1(t,y)
y
dy = (t*D*y+A*t+C+sqrt(B*y))/((t*t-t)*D)
end
end
Paul
2024-1-7
Yes, I thought we already established that original code eventually leads to an integration step where y becomes negative and dy becomes complex valued, at which point y will have a non-zero imaginary component. So I'm still not sure what the concern is with the original code.
Maybe you're just referring to the display of the output in the Matlab window, which depends on the format of the output display and how Matlab diplays 0 values (for either real or imaginary parts) of complex-valued arrays.
Consider the real array
x = [0;eps]
x = 2×1
1.0e-15 *
0
0.2220
Notice that with the defaut formating the first element is displayed as "0", which indicates that it's exactly zero ( and multiplied by 1e-15)
But, with a complex array
y = [0; 0 + 1i*eps; eps]
y =
1.0e-15 *
0.0000 + 0.0000i
0.0000 + 0.2220i
0.2220 + 0.0000i
Ther real and imaginary parts of elements where those parts are exactly zero are displayed with "0.0000".
But if we change format, we see those parts are, in fact, exactly zero
format long e
y
y =
0.000000000000000e+00 + 0.000000000000000e+00i
0.000000000000000e+00 + 2.220446049250313e-16i
2.220446049250313e-16 + 0.000000000000000e+00i
Sam Chak
2024-1-7
编辑：Sam Chak
2024-1-7
Hi @Kaiwei
As noted by other users, the complex-valued solution is caused by the square root term 'sqrt(B*y)' when the output y is negative. In theory, the integration will yield a negative output y when y-prime (
) is negative. However, the output y vector in the numerical solution does not have a negative real part. Thus, in a human sense, the argument in the square root term should remain positive, as you would expect a monotonically increasing pattern in the output.
) is negative. However, the output y vector in the numerical solution does not have a negative real part. Thus, in a human sense, the argument in the square root term should remain positive, as you would expect a monotonically increasing pattern in the output.Here, using the conventional Runge-Kutta 4th-order (RK4) solver, we can examine what really happens to the initial step when the step size is very small. The solution of RK4 consists of 4 terms, 
, 
, 
, and 
. The result below shows that 
becomes complex-valued because become negative.
, , , and . The result below shows that becomes complex-valued because become negative.format long g
%% Kaiwei's model
% Parameters
A = -4.32;
B = 4;
C = 1.2;
D = 0.8;
a = 5/18; % 0.277777777777777777777777777777777
% Ordinary Differential Equations
odefcn = @(t, y) (t*D*y + A*t + C + sqrt(B*y))/(D*(t.^2 - t));
% Simulation time span
tStart = a; % 5/18
h = 0.00000013/18; % a very small step size
tEnd = a + h; % check what happen to the initial step
% tEnd = 1 - 0.0013/18; % singularity occurs at 18/18 = 1
t = tStart:h:tEnd;
% Initial value
y0 = 0;
% Call RK4 Solver
yRK4 = RK4Solver(odefcn, t, y0);
k1 =
1.38350869222519e-15
k2 =
9.71721452549591e-08
k3 =
-1.36233392341977e-07
k4 =
1.94399999139241e-07 - 3.90884873194043e-07i
y =
0 + 0i 1.39963664796411e-16 - 4.70509568735031e-16i
% Plotting the solutions
plot(t, yRK4, '-o')
Warning: Imaginary parts of complex X and/or Y arguments ignored.
grid on
xlabel('\alpha')
ylabel('U')
%% Runge-Kutta 4th-order Solver
function y = RK4Solver(f, x, y0)
y(:, 1) = y0; % initial condition
h = x(2) - x(1); % step size
n = length(x); % number of steps
for j = 1 : n-1
k1 = f(x(j), y(:, j))
k2 = f(x(j) + h/2, y(:, j) + h/2*k1)
k3 = f(x(j) + h/2, y(:, j) + h/2*k2)
k4 = f(x(j) + h, y(:, j) + h*k3)
y(:, j+1) = y(:, j) + h/6.0*(k1 + 2*k2 + 2*k3 + k4)
end
end
3 个评论
Sam Chak
2024-1-7
Hi @Kaiwei
An implicit ODE doesn't look like yours. Essentially, as @Torsten explained, once the complex-valued solution appears in y, the imaginary part 'cannot be removed'. The values of the imaginary part do not contribute significantly because they are very, very small. You can evaluate the values of the imaginary part below:
By the way, just out of curiosity, what physical phenomenon does your differential equation describe?
format long g
%% Kaiwei's model
% Parameters
A = -4.32;
B = 4;
C = 1.2;
D = 0.8;
a = 5/18; % 0.277777777777777777777777777777777
% Ordinary Differential Equations
odefun = @(t, y) (t*D*y + A*t + C + sqrt(B*y))/(D*(t.^2 - t));
% Simulation time span
tStart = a; % 5/18
h = 0.13/18; % a very small step size
tEnd = 1 - 0.13/18; % singularity occurs at 18/18 = 1
t = tStart:h:tEnd;
% Call ode45 solver
tspan = a:0.13/18:(1 - 0.13/18);
y0 = 0;
options = odeset('Reltol', 1e-10, 'MaxStep', (13/18)/1e3);
sol = ode45(odefun, tspan, y0, options);
[y, yp] = deval(sol, tspan) % yp is y-prime (y')
y =
Columns 1 through 3
0 + 0i 0.000149371117698822 - 8.74290466347223e-11i 0.000593850411598535 - 7.10875020470179e-12i
Columns 4 through 6
0.00132819985638987 - 1.65895186942128e-12i 0.00234745375139398 - 5.95623798827648e-13i 0.00364690051279341 - 2.70617263766076e-13i
Columns 7 through 9
0.00522206707631183 - 1.42626355904756e-13i 0.00706870451694508 - 8.32449774310712e-14i 0.00918277488975191 - 5.23372549180808e-14i
Columns 10 through 12
0.0115604391731636 - 3.48185724770688e-14i 0.0141980462065993 - 2.42143209444681e-14i 0.017092122526636 - 1.7451314136623e-14i
Columns 13 through 15
0.0202393630170698 - 1.2950899838263e-14i 0.0236366222978676 - 9.84872225180023e-15i 0.0272809067864375 - 7.64594808207701e-15i
Columns 16 through 18
0.031169367372005 - 6.04168931127621e-15i 0.0352992926503299 - 4.84749674404772e-15i 0.0396681026716628 - 3.94145189756681e-15i
Columns 19 through 21
0.0442733431598302 - 3.24241862736225e-15i 0.0491126801647401 - 2.69504368690119e-15i 0.0541838951145055 - 2.26072375099067e-15i
Columns 22 through 24
0.0594848802368446 - 1.91200343568839e-15i 0.0650136343225046 - 1.62900906703916e-15i 0.0707682588062083 - 1.39712457246199e-15i
Columns 25 through 27
0.076746954143092 - 1.20544451978269e-15i 0.0829480164608167 - 1.04572453300378e-15i 0.089369834469532 - 9.11656647100207e-16i
Columns 28 through 30
0.0960108866136756 - 7.98360967212335e-16i 0.10286973845123 - 7.02023811354232e-16i 0.109945040247545 - 6.19636631824521e-16i
Columns 31 through 33
0.117235524772198 - 5.4880528800622e-16i 0.124740005288621 - 4.87609096054951e-16i 0.132457373727367 - 4.34495541141578e-16i
Columns 34 through 36
0.140386599034969 - 3.88200839026214e-16i 0.148526725691339 - 3.4768943846025e-16i 0.156876872389597 - 3.12107543817558e-16i
Columns 37 through 39
0.165436230873105 - 2.80747114809518e-16i 0.174204064925331 - 2.53017765759942e-16i 0.183179709508972 - 2.28424671261845e-16i
Columns 40 through 42
0.192362570051574 - 2.06551075043324e-16i 0.201752121875624 - 1.87044353132282e-16i 0.2113479097719 - 1.69604840881079e-16i
Columns 43 through 45
0.22114954771557 - 1.53976823628015e-16i 0.231156718725351 - 1.39941231905414e-16i 0.241369174866789 - 1.27309687633596e-16i
Columns 46 through 48
0.251786737401525 - 1.15919627230175e-16i 0.262409297085286 - 1.05630287860364e-16i 0.27323681461817 - 9.6319389096335e-17i
Columns 49 through 51
0.284269321251758 - 8.78803776355824e-17i 0.295506919558587 - 8.02201300829646e-17i 0.306949784370563 - 7.3257030072572e-17i
Columns 52 through 54
0.318598163894094 - 6.69193526381729e-17i 0.330452381010976 - 6.11439018154854e-17i 0.34251283477549 - 5.58748577892358e-17i
Columns 55 through 57
0.35478000211971 - 5.10627980989256e-17i 0.367254439780776 - 4.66638639585549e-17i 0.379936786465816 - 4.26390479865406e-17i
Columns 58 through 60
0.392827765272418 - 3.89535838593367e-17i 0.405928186384993 - 3.55764218098556e-17i 0.41923895007025 - 3.24797766561675e-17i
Columns 61 through 63
0.432761049998147 - 2.96387372971496e-17i 0.446495576918451 - 2.7030928451849e-17i 0.460443722727234 - 2.46362169287613e-17i
Columns 64 through 66
0.474606784962578 - 2.24364559537015e-17i 0.488986171774456 - 2.04152621110871e-17i 0.503583407420437 - 1.85578203036515e-17i
Columns 67 through 69
0.518400138346644 - 1.68507128422229e-17i 0.533438139922607 - 1.52817693662587e-17i 0.548699323909487 - 1.38399347882945e-17i
Columns 70 through 72
0.564185746754093 - 1.25151528683543e-17i 0.579899618816521 - 1.12982633714691e-17i 0.595843314657762 - 1.01809110540147e-17i
Columns 73 through 75
0.61201938453603 - 9.15546497179623e-18i 0.62843056728772 - 8.21494681221603e-18i 0.645079804802189 - 7.35296713065195e-18i
Columns 76 through 78
0.661970258340474 - 6.56366852246032e-18i 0.679105326998836 - 5.84167489100917e-18i 0.696488668681494 - 5.18204608234435e-18i
Columns 79 through 81
0.714124224026862 - 4.58023725138386e-18i 0.73201624383333 - 4.03206240532993e-18i 0.750169320661215 - 3.53366163927338e-18i
Columns 82 through 84
0.768588425456956 - 3.08147163839017e-18i 0.787278950268152 - 2.67219907204337e-18i 0.80624675841374 - 2.30279654859767e-18i
Columns 85 through 87
0.825498243872288 - 1.97044083667637e-18i 0.845040402197028 - 1.67251308955588e-18i 0.864880916025902 - 1.40658083475474e-18i
Columns 88 through 90
0.885028259332959 - 1.1703815106714e-18i 0.905491826130283 - 9.61807345959525e-19i 0.926282091651604 - 7.7889138413146e-19i
Columns 91 through 93
0.947410817595949 - 6.19794453452836e-19i 0.968891318607491 - 4.82792866300509e-19i 0.990738816343995 - 3.66266594645507e-19i
Columns 94 through 96
1.01297092322675 - 2.68687592957645e-19i 1.03560832649641 - 1.88607792480594e-19i 1.05867579859296 - 1.24645994133376e-19i
Columns 97 through 99
1.0822037773592 - 7.5472247421798e-20i 1.10623104017854 - 3.97867636832611e-20i 1.13080979181703 - 1.62860379486857e-20i
Column 100
1.15601747352413 - 3.59242447518241e-21i
yp =
Columns 1 through 3
1.38350869222519e-15 + 0i 0.0412369899175437 + 4.40005118066934e-08i 0.0817288630656546 + 1.77304473973381e-09i
Columns 4 through 6
0.121514896524142 + 2.73607498895247e-10i 0.160631338173492 + 7.31318522825841e-11i 0.199112267007545 + 2.6403988112413e-11i
Columns 7 through 9
0.236989683230474 + 1.15271713733352e-11i 0.274293683577695 + 5.73608902422251e-12i 0.311052621800149 + 3.14080850813915e-12i
Columns 10 through 12
0.347293253472141 + 1.84981630917606e-12i 0.383040866745612 + 1.15382448742241e-12i 0.418319400659591 + 7.5382585698402e-13i
Columns 13 through 15
0.453151552430337 + 5.11654572928551e-13i 0.4875588749681 + 3.58563302684458e-13i 0.52156186570885 + 2.58194406800592e-13i
Columns 16 through 18
0.555180047713462 + 1.90310550568753e-13i 0.588432043869964 + 1.43145856790338e-13i 0.621335644934 + 1.09598088648478e-13i
Columns 19 through 21
0.653907872055901 + 8.5237673679364e-14i 0.686165034367851 + 6.72214303215872e-14i 0.718122782139676 + 5.36773865080462e-14i
Columns 22 through 24
0.749796155955552 + 4.33447220540839e-14i 0.781199632314842 + 3.53567615502601e-14i 0.812347166017644 + 2.91067024652601e-14i
Columns 25 through 27
0.843252229658484 + 2.41625445856791e-14i 0.873927850518993 + 2.02120231724867e-14i 0.904386645122159 + 1.70262608107417e-14i
Columns 28 through 30
0.934640851685789 + 1.44353427646046e-14i 0.964702360691227 + 1.23116333428449e-14i 0.994582743764448 + 1.05582065771976e-14i
Columns 31 through 33
1.02429328105009 + 9.10071075663193e-15i 1.05384498724463 + 7.88157318792078e-15i 1.08324863644258 + 6.85582209559713e-15i
Columns 34 through 36
1.11251478593874 + 5.98804048584521e-15i 1.14165379912041 + 5.25012194189341e-15i 1.17067586757609 + 4.61960096204647e-15i
Columns 37 through 39
1.19959103254012 + 4.07839928533953e-15i 1.22840920578838 + 3.61187639861934e-15i 1.25714019009539 + 3.20810454805129e-15i
Columns 40 through 42
1.28579369936101 + 2.85731090770287e-15i 1.31437937851245 + 2.55144524111499e-15i 1.34290682328711 + 2.28384251692962e-15i
Columns 43 through 45
1.37138560000185 + 2.04895790591941e-15i 1.39982526541584 + 1.84215734205711e-15i 1.42823538679662 + 1.65955102342642e-15i
Columns 46 through 48
1.45662556230317 + 1.49786030837365e-15i 1.48500544180443 + 1.35431074134438e-15i 1.51338474825848 + 1.22654564171525e-15i
Columns 49 through 51
1.5417732997859 + 1.11255596409401e-15i 1.57018103258012 + 1.01062310202328e-15i 1.59861802480994 + 9.19272039575295e-16i
Columns 52 through 54
1.62709452168271 + 8.37232815670108e-16i 1.65562096185387 + 7.63408697048546e-16i 1.68420800538754 + 6.9684978932963e-16i
Columns 55 through 57
1.71286656349555 + 6.36731074946768e-16i 1.74160783030945 + 5.82334069494163e-16i 1.77044331697105 + 5.33031447263108e-16i
Columns 58 through 60
1.79938488836392 + 4.88274112424406e-16i 1.82844480285163 + 4.47580291943331e-16i 1.85763575543908 + 4.10526305639917e-16i
Columns 61 through 63
1.88697092483386 + 3.7673873223212e-16i 1.91646402495494 + 3.45887741115285e-16i 1.946129361521 + 3.17681400664937e-16i
Columns 64 through 66
1.97598189445068 + 2.91860807043591e-16i 2.0060373069274 + 2.68195904441786e-16i 2.03631208212506 + 2.46481889640262e-16i
Columns 67 through 69
2.06682358876532 + 2.26536111730079e-16i 2.09759017688673 + 2.08195392550054e-16i 2.12863128546264 + 1.91313705514866e-16i
Columns 70 through 72
2.15996756381745 + 1.75760160505783e-16i 2.19162100917675 + 1.61417250773342e-16i 2.22361512316354 + 1.48179324673976e-16i
Columns 73 through 75
2.25597509064679 + 1.35951250785177e-16i 2.2887279850948 + 1.24647249722392e-16i 2.32190300552654 + 1.14189869982124e-16i
Columns 76 through 78
2.35553175135414 + 1.04509088495717e-16i 2.3896485429493 + 9.55415194077157e-17i 2.42429079776133 + 8.72297169825614e-17i
Columns 79 through 81
2.45949947442504 + 7.95215605686087e-17i 2.49531960074745 + 7.23697112710573e-17i 2.53180090607608 + 6.57311314574805e-17i
Columns 82 through 84
2.56899858479841 + 5.95666594851577e-17i 2.60697422629152 + 5.38406331366313e-17i 2.64579695857514 + 4.85205562131107e-17i
Columns 85 through 87
2.68554486980932 + 4.35768035971041e-17i 2.72630679611321 + 3.89823608896301e-17i 2.76818459995312 + 3.47125954922032e-17i
Columns 88 through 90
2.81129611712182 + 3.07450567897552e-17i 2.85577903323742 + 2.7059303970948e-17i 2.90179608222685 + 2.36367611145594e-17i
Columns 91 through 93
2.94954217495992 + 2.04606006792668e-17i 2.99925443371145 + 1.75156588395919e-17i 3.05122676315099 + 1.47883899673817e-17i
Columns 94 through 96
3.1058318206262 + 1.22668745670366e-17i 3.16355572332299 + 9.94090880566564e-18i 3.22505622698764 + 7.8022336883866e-18i
Columns 97 through 99
3.29126820348546 + 5.8450341876331e-18i 3.36361695552618 + 4.06704084083267e-18i 3.44452726475602 + 2.47226704657304e-18i
Column 100
3.53905122463319 + 1.07990810763503e-18i
plot(tspan, y, tspan, yp), grid
Warning: Imaginary parts of complex X and/or Y arguments ignored.
xlabel('\alpha'), legend('y', 'y-prime')
Kaiwei
2024-1-7
Seems I can not solve this problem? Do you think that will be a big problem?
Actually I study a bandit learning model (economics theory field), and this U stands for my value function. The belief updates in a Bayesian way.
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!发生错误
由于页面发生更改，无法完成操作。请重新加载页面以查看其更新后的状态。
Translated by 
您也可以从以下列表中选择网站：
如何获得最佳网站性能
选择中国网站（中文或英文）以获得最佳网站性能。其他 MathWorks 国家/地区网站并未针对您所在位置的访问进行优化。
美洲
- América Latina (Español)
- Canada (English)
- United States (English)
欧洲
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom(English)
亚太
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)