Getting complex solution with Matlab ode45

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);

[Paul]

What is the desired solution?

[Torsten]

t^2-t is negative for 0 < t < 1. Since you start with an initial value for y not far from 0, it's not surprising that dy < 0 leads to y < 0 and thus a complex square root in the course of the integration.

[Kaiwei]

The desired solution should be a monotonically increasing one starts from a, and smoothly intersects with horizontal line at a.

[Torsten]

This will not happen - the denominator is negative for 0 < t < 1 and thus your solution will always be decreasing for your choice of parameters.

[Kaiwei]

Actually at around a, the numerator is also negative since A*a+C is negative. So around a, y' should be positive.

[Paul]

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]

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]

About the monotonicity, I get that from maximum principle so I guess it's for sure.

[Torsten]

So you deny the facts ? @Paul's call to "odefun_mon_pc1" clearly shows that your function is decreasing, and as soon as y becomes negative, sqrt(B*y) becomes complex-valued. If a>1, it's another story.

[Kaiwei]

I actually want y(a) to be 0, since I feel it's possible to cause problems in sqrt, I use 0.0001 instead.

Answer 1

[Paul]

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.

[Kaiwei]

By the way how can I get a callable function y and y'?

(for example use sol=ode45() and define y=@(t) deval(sol,t) but it's not defined on everywhere, so sometimes it's not convenient)

[Kaiwei]

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]

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]

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]

@Kaiwei

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]

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]

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?

[Kaiwei]

Thanks for the second part. That's very helpful

[Torsten]

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]

@Kaiwei

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]

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]

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]

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]

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]

@Kaiwei,

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

Answer 2

[Sam Chak]

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.

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.

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

[Kaiwei]

can you say something more about this comment? I didn't get how to solve an implicit ODE.

[Sam Chak]

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]

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.