denitsastaicova

50 Reputation

6 Badges

12 years, 153 days

MaplePrimes Activity


These are replies submitted by denitsastaicova

Ok, here's an extremely simplified example of my system which I rounded to 3 digits of precision.

A very important difference with the non-rounded system is that in it, the imaginary part of SQ (and also of x(t)) is 10^(-17) and not 10^(-7). But otherwise, you see exactly the same problem. If you check out the imaginary part of the term I think the problem is - SQ, before the so-called "singularity" point it's 10^(-6), and after, it's 10^(-5). So it's not like we see some huge jump in the imaginary part. In the non-rounded system it's 10^(-18) and 10^(-12) respectively. In the real system (non-rounded, non-simplified), it's changes between 10^(-17)*I and 0*I, which looks very safe to me.
 

restart:

Digits:=15:

Eq:=3.000*x(t)^2*diff(y(t), t)*diff(x(t), t) + 1.000*x(t)^3*diff(y(t), t, t) = 0, diff(x(t), t) = sqrt(6)*sqrt(0.500*diff(y(t), t)^2 + ((5.000/(x(t)^3) + sqrt(-37.037 + 25.000/(x(t)^6)))^(1/3) + 3.333/(5.000/(x(t)^3) + sqrt(-37.037 + 25.000/(x(t)^6)))^(1/3))^2/4 + 3*((5.000/(x(t)^3) + sqrt(-37.037 + 25.000/(x(t)^6)))^(1/3) + 3.333/(5.000/(x(t)^3) + sqrt(-37.037 + 25.000/(x(t)^6)))^(1/3))/(4*x(t)^3))*x(t)/6

3.000*x(t)^2*(diff(y(t), t))*(diff(x(t), t))+1.000*x(t)^3*(diff(diff(y(t), t), t)) = 0, diff(x(t), t) = (1/6)*6^(1/2)*(.500*(diff(y(t), t))^2+(1/4)*((5.000/x(t)^3+(-37.037+25.000/x(t)^6)^(1/2))^(1/3)+3.333/(5.000/x(t)^3+(-37.037+25.000/x(t)^6)^(1/2))^(1/3))^2+(3/4)*((5.000/x(t)^3+(-37.037+25.000/x(t)^6)^(1/2))^(1/3)+3.333/(5.000/x(t)^3+(-37.037+25.000/x(t)^6)^(1/2))^(1/3))/x(t)^3)^(1/2)*x(t)

(1)

IC:=y(0) = 1, D(y)(0) = 0, x(0) = 1/10

y(0) = 1, (D(y))(0) = 0, x(0) = 1/10

(2)

 

f_int:=1.5:
Ti:=time():sol1:=dsolve(evalf({Eq,IC}), numeric,complex=true): A1:=sol1(1);A2:=sol1(f_int);

[t = 1., x(t) = 1.31705692386913-0.100091386113909e-5*I, y(t) = 1.00000000000000+0.*I, diff(y(t), t) = 0.+0.*I]

 

[t = 1.5, x(t) = 1.91197228324010-0.122178940697193e-5*I, y(t) = 1.00000000000000+0.*I, diff(y(t), t) = 0.+0.*I]

(3)

sol1(0.625)

[t = .625, x(t) = .936467027370079+0.905489928187092e-7*I, y(t) = .999999999999996+0.*I, diff(y(t), t) = 0.+0.*I]

(4)

sol1(0.627)

[t = .627, x(t) = .938441466071138+0.856451296997087e-7*I, y(t) = 1.00000000000000+0.*I, diff(y(t), t) = 0.+0.*I]

(5)

SQ:=((5.000/(x(t)^3) + sqrt(-37.037 + 25.000/(x(t)^6)))^(1/3) + 3.333/(5.000/(x(t)^3) + sqrt(-37.037 + 25.000/(x(t)^6)))^(1/3))^2

((5.000/x(t)^3+(-37.037+25.000/x(t)^6)^(1/2))^(1/3)+3.333/(5.000/x(t)^3+(-37.037+25.000/x(t)^6)^(1/2))^(1/3))^2

(6)

subs((4),SQ)

13.3332040917100-0.864213791841415e-6*I

(7)

subs((5),SQ)

13.3145138736522-0.489002159265053e-4*I

(8)

f_int:=1.5:
Ti:=time():sol1:=dsolve(evalf({Eq,IC}), numeric): A1:=sol1(1);A2:=sol1(f_int);

Error, (in sol1) cannot evaluate the solution further right of .62512724, probably a singularity

 

Error, (in sol1) cannot evaluate the solution further right of .62512724, probably a singularity

 

 

 


 

Download example_ode.mw

@mmcdara 

Hi, thanks for replying. Unfortunately, trying to separate the real and imaginary part so far leads only to a huge memory leak and well, had to be killed. So yeah, that seems not to work. The problem is that my system consists of 2 coupled ODEs and because of the sqrt, it is an implicit one. I think I see where the imaginary part comes from (solving a cubic equation), but any attempt to fix it so far hasn't worked. Adding "Re" to the parameters I think might be problematic leads to

"Could not convert to an explicit first order system due to 'RootOf'" or random complains about Re or abs.

 Still, thanks for the example! If anyone can think of some clever way to add Re(f(x)) to an ODE, which won't crash the solver, I'd be very grateful.

Yeah, just %f should be %e:

printf("%.5Ze %.5e %g",-0.123456e-18-0.123456e-15*I, 0.1334423423*10^(-15),3);

Otherwise the %f will cut the floating parts of the comples number.

Thank you, Alejandro! You and pagan really helped me :)

So, for the record, the overall routine should be something like:

A:=[-0.9123456e-9-0.123456e-17*I,0.1334423423*10^(-15),3]:
fopen(fi1,APPEND);
fprintf(fi1,"%.5Zg %.5e %g\n",op(A)):
fclose(fi1);

I guess my question was kind of trivial for the people who has programmed on C, but for me, it was a nightmare to do it without help. So, thanks.

Yeah, just %f should be %e:

printf("%.5Ze %.5e %g",-0.123456e-18-0.123456e-15*I, 0.1334423423*10^(-15),3);

Otherwise the %f will cut the floating parts of the comples number.

Thank you, Alejandro! You and pagan really helped me :)

So, for the record, the overall routine should be something like:

A:=[-0.9123456e-9-0.123456e-17*I,0.1334423423*10^(-15),3]:
fopen(fi1,APPEND);
fprintf(fi1,"%.5Zg %.5e %g\n",op(A)):
fclose(fi1);

I guess my question was kind of trivial for the people who has programmed on C, but for me, it was a nightmare to do it without help. So, thanks.

By the way, any idea how to make the 3 numbers to be written in the same line and then to break to a new line.

The second part is easier since fprintf(f,`%a`,x,`%\n`);  in the integer part of the procedure will do the work. But how can I get the other two on the same line? I couldn't find in the help how to avoid the new line when fprintf is invoked on the top level.
 

By the way, any idea how to make the 3 numbers to be written in the same line and then to break to a new line.

The second part is easier since fprintf(f,`%a`,x,`%\n`);  in the integer part of the procedure will do the work. But how can I get the other two on the same line? I couldn't find in the help how to avoid the new line when fprintf is invoked on the top level.
 

Yes, this works very nicely for the example I provided, however for:

A:=[-0.9123456e-9-0.123456e-17*I,0.1334423423*10^(-15),3]: writedata('terminal',A,string,T);
-0.00000-0.00000I
1.3344e-16
3

So with little experimenting, I found that

T:=proc(f,x) if type(x,float) then fprintf(f,`%.5g`,x); elif type(x,complex(float)) then fprintf(f,`%.5Ze`,x); else fprintf(f,`%a`,x); end if; NULL; end proc:

and

 T:=proc(f,x) if type(x,float) then fprintf(f,`%.5g`,x); elif type(x,complex(float)) then fprintf(f,`%.5Zg`,x); else fprintf(f,`%a`,x); end if; NULL; end proc:

would print precisely the desired:

-9.12346e-10-1.23456e-18I
1.3344e-16
3
 

Thanks a lot for the help!

Yes, this works very nicely for the example I provided, however for:

A:=[-0.9123456e-9-0.123456e-17*I,0.1334423423*10^(-15),3]: writedata('terminal',A,string,T);
-0.00000-0.00000I
1.3344e-16
3

So with little experimenting, I found that

T:=proc(f,x) if type(x,float) then fprintf(f,`%.5g`,x); elif type(x,complex(float)) then fprintf(f,`%.5Ze`,x); else fprintf(f,`%a`,x); end if; NULL; end proc:

and

 T:=proc(f,x) if type(x,float) then fprintf(f,`%.5g`,x); elif type(x,complex(float)) then fprintf(f,`%.5Zg`,x); else fprintf(f,`%a`,x); end if; NULL; end proc:

would print precisely the desired:

-9.12346e-10-1.23456e-18I
1.3344e-16
3
 

Thanks a lot for the help!

Closer, but as you notice the second number is still far from the 5 digit goal. I really don't understand why it's so hard to do this. Maybe I miss the C programming knowledge, but still...

But your variant is at least fixing the problem for the complex number.

Thanks!

Closer, but as you notice the second number is still far from the 5 digit goal. I really don't understand why it's so hard to do this. Maybe I miss the C programming knowledge, but still...

But your variant is at least fixing the problem for the complex number.

Thanks!

Yup, I want that :) Especially since I prepare the files on one pc and then run them on another. Thanks for letting me know of this option.

Yeah, the power of procrastination :)

I love this site! It reminds me very much of my life right now :)

I want something like: p:=1; for i from 1 by p to 100 do A:=evalf(F(p)); if A <1 p:=p/2;fi; if A >1 p:=p*2 fi; end do; The idea is to modify the step according to my needs. But if I make: p:=1; for i from 1 by p to 100 do A:=evalf(F(p)); if A<1 p:=p/2;fi; if A>1 p:=p*2 fi; print(A,i,p) end do; I see that p changes as required, but i doesn't-I see something like: 0.5, 1, 1, 1, 2, 2 1,3, 4 I'm just creating the example now, but the idea is that p changes with the if clause, while the step in the for-statement which is supposed do be p isn't changing-it's always 1-what was set in the beginning. Anyway, I found a way around this setting a different counter for the for-cycle and using while. It works for some simply examples, I hope it will work with the more complex ones too.
I want something like: p:=1; for i from 1 by p to 100 do A:=evalf(F(p)); if A <1 p:=p/2;fi; if A >1 p:=p*2 fi; end do; The idea is to modify the step according to my needs. But if I make: p:=1; for i from 1 by p to 100 do A:=evalf(F(p)); if A<1 p:=p/2;fi; if A>1 p:=p*2 fi; print(A,i,p) end do; I see that p changes as required, but i doesn't-I see something like: 0.5, 1, 1, 1, 2, 2 1,3, 4 I'm just creating the example now, but the idea is that p changes with the if clause, while the step in the for-statement which is supposed do be p isn't changing-it's always 1-what was set in the beginning. Anyway, I found a way around this setting a different counter for the for-cycle and using while. It works for some simply examples, I hope it will work with the more complex ones too.
hi, i have similar problem in linux version of Maple 11. In order to make it start at all, i had to rename 2 folders, as i'm on 64bit AMD. And now, Maple starts but says - no kernel connection. I'm pretty sure it's not the firewall, so i have no idea what to do. Tried everything i could come up with. And ideas?
1 2 Page 1 of 2