rtoy / maxima

A Clone of Maxima's repo
Other
0 stars 0 forks source link

cd(u, m) vs sn(u + K(m), m) #2263

Open rtoy opened 3 months ago

rtoy commented 3 months ago

Imported from SourceForge on 2024-07-05 22:36:19 Created by vladutzanu on 2016-09-14 14:52:52 Original: https://sourceforge.net/p/maxima/bugs/3214


Hello

Given the following input:

q(x, y) := exp(-%pi * elliptic_kc( ( 1 - x^y )^(1/y) ) / elliptic_kc(x))$
Theta[2](x, y, n) := 2 * y^0.25 * sum( y^(k * (k + 1)) * cos((2 * k + 1) * x), k, 0, n );
Theta[3](x, y, n) := 1 + 2 * sum( y^(k^2) * cos(2 * k * x), k, 1, n );
[_u, _m, N] : [25, 0.999999, 5]$
U(x, y) := %pi * x / (2 * elliptic_kc(y))$
M(y) := q(y, 1)$
float(Theta[2](U(_u, _m), M(_m), N) / Theta[3](U(_u, _m), M(_m), N) / _m^0.25);
jacobi_cd(_u, _m);
plot2d([Theta[2](U(x, _m), M(_m), N) / Theta[3](U(x, _m), M(_m), N) / _m^0.25 - jacobi_cd(x + 0*elliptic_kc(_m), _m)],
    [x, 0, 4*elliptic_kc(_m)])$

I ran into some strange results when trying to approximate jacobi_cd() with theta functions. The errors are hinted as early as _m=0.9, and get bigger from there on. Using jacobi_sn( u + elliptic_kc(m), m ) the results are correct. It's true it's extremely unlikely (I'd say impossible) to get to such a modulus in my quests, but I thought I'd let you know.

rtoy commented 3 months ago

Imported from SourceForge on 2024-07-05 22:36:20 Created by robert_dodier on 2016-09-16 05:24:28 Original: https://sourceforge.net/p/maxima/bugs/3214/#1f36


Diff:


--- old
+++ new
@@ -2,15 +2,17 @@

 Given the following input:

-> q(x, y) := exp(-%pi * elliptic_kc( ( 1 - x^y )^(1/y) ) / elliptic_kc(x))$
-> Theta[2](x, y, n) := 2 * y^0.25 * sum( y^(k * (k + 1)) * cos((2 * k + 1) * x), k, 0, n );
-> Theta[3](x, y, n) := 1 + 2 * sum( y^(k^2) * cos(2 * k * x), k, 1, n );
-> [_u, _m, N] : [25, 0.999999, 5]$
-> U(x, y) := %pi * x / (2 * elliptic_kc(y))$
-> M(y) := q(y, 1)$
-> float(Theta[2](U(_u, _m), M(_m), N) / Theta[3](U(_u, _m), M(_m), N) / _m^0.25);
-> jacobi_cd(_u, _m);
-> wxplot2d([Theta[2](U(x, _m), M(_m), N) / Theta[3](U(x, _m), M(_m), N) / _m^0.25 - jacobi_cd(x + 0*elliptic_kc(_m), _m)],
->     [x, 0, 4*elliptic_kc(_m)])$
+~~~~
+q(x, y) := exp(-%pi * elliptic_kc( ( 1 - x^y )^(1/y) ) / elliptic_kc(x))$
+Theta[2](x, y, n) := 2 * y^0.25 * sum( y^(k * (k + 1)) * cos((2 * k + 1) * x), k, 0, n );
+Theta[3](x, y, n) := 1 + 2 * sum( y^(k^2) * cos(2 * k * x), k, 1, n );
+[_u, _m, N] : [25, 0.999999, 5]$
+U(x, y) := %pi * x / (2 * elliptic_kc(y))$
+M(y) := q(y, 1)$
+float(Theta[2](U(_u, _m), M(_m), N) / Theta[3](U(_u, _m), M(_m), N) / _m^0.25);
+jacobi_cd(_u, _m);
+plot2d([Theta[2](U(x, _m), M(_m), N) / Theta[3](U(x, _m), M(_m), N) / _m^0.25 - jacobi_cd(x + 0*elliptic_kc(_m), _m)],
+    [x, 0, 4*elliptic_kc(_m)])$
+~~~~

-I ran into some strange results when trying to approximate jacobi_cd() with theta functions. The errors are hinted as early as _m=0.9, and get bigger from there on. Using jacobi_sn( u + elliptic_kc(m), m ) the results are correct. It's true it's extremely unlikely (I'd say impossible) to get to such a modulus in my quests, but I thought I'd let you know.
+I ran into some strange results when trying to approximate jacobi_cd() with theta functions. The errors are hinted as early as `_m=0.9`, and get bigger from there on. Using jacobi_sn( u + elliptic_kc(m), m ) the results are correct. It's true it's extremely unlikely (I'd say impossible) to get to such a modulus in my quests, but I thought I'd let you know.
rtoy commented 3 months ago

Imported from SourceForge on 2024-07-05 22:36:24 Created by robert_dodier on 2016-09-16 05:24:29 Original: https://sourceforge.net/p/maxima/bugs/3214/#eb70


Put in formatting for code. Also changed wxplot2d to plot2d so that the example works for command line Maxima.

rtoy commented 3 months ago

Imported from SourceForge on 2024-07-05 22:36:27 Created by vladutzanu on 2016-09-17 14:17:47 Original: https://sourceforge.net/p/maxima/bugs/3214/#3b73


I thought I'd check and sure enough, I forgot to specify that the q(x,y) is the nome q with squared modulus x (y=1), where x'=1-x, or normal modulus (y=2), where x'=sqrt(1-x^2). Sorry for the omission.

rtoy commented 3 months ago

Imported from SourceForge on 2024-07-05 22:36:31 Created by robert_dodier on 2016-09-17 20:36:30 Original: https://sourceforge.net/p/maxima/bugs/3214/#3b73/9d49


Do you mean to say that it's not a bug after all? If so, can you please close the report. If indeed it is still a bug, can you please update the description to include the new definition so that someone can just copy and paste the example and get the same result.

rtoy commented 3 months ago

Imported from SourceForge on 2024-07-05 22:36:34 Created by vladutzanu on 2016-09-20 19:57:23 Original: https://sourceforge.net/p/maxima/bugs/3214/#3b73/9d49/eb72


I'm sorry for the delay. No, that was just a specification for that particular implementation of the nome q. The bug is according to the discussion with the title "cd(u,m) vs sn(u+K(m),m)", found in the maxima-discuss mailing list.

rtoy commented 3 months ago

Imported from SourceForge on 2024-07-05 22:36:38 Created by robert_dodier on 2016-09-27 15:25:57 Original: https://sourceforge.net/p/maxima/bugs/3214/#3754


OK, I ran the code in the description and here is the plot that I get. Do I understand correctly that the problem is that the plotted curve should be exactly zero? If so, I guess the problem is that the numerical approximation is not accurate throughout the range of x.

Attachments:

rtoy commented 3 months ago

Imported from SourceForge on 2024-07-05 22:36:42 Created by vladutzanu on 2016-09-27 16:34:07 Original: https://sourceforge.net/p/maxima/bugs/3214/#3754/d4db


Yes, that's the problem. Using jacobi_sn(u+elliptic_kc(m),m) seems to give correct results, i.e. I can only see the errors from the approximation with the theta functions.

rtoy commented 3 months ago

Imported from SourceForge on 2024-07-05 22:36:45 Created by rtoy on 2016-09-29 05:02:17 Original: https://sourceforge.net/p/maxima/bugs/3214/#2447


Not sure what the problem is. I tried plot2d(jacobi_cd(u,.9)-jacobi_sn(u+elliptic_kc(.9),.9),[u,0,4*elliptic_kc(.9)]); and the plot shows a max diff of about 3.5e-14. It seems like they're identical to within rounding errors.

rtoy commented 3 months ago

Imported from SourceForge on 2024-07-05 22:36:49 Created by vladutzanu on 2016-09-29 05:43:24 Original: https://sourceforge.net/p/maxima/bugs/3214/#2447/8798


The modulus is not small enough. Add one 9 to make it 0.99 and the errors will increase. Add some more, and you'll see what I mean. For my application there, I would need at least 0.99997671.

rtoy commented 3 months ago

Imported from SourceForge on 2024-07-05 22:36:52 Created by rtoy on 2016-10-04 02:50:58 Original: https://sourceforge.net/p/maxima/bugs/3214/#b576


Ah, ok. 0.99998 really shows a difference of about 1e-6, which is rather large. As best as I can tell, the problem is that jacobi_cn is not very accurate for large modulus and argument. I need to find a different way of computing it.