rtoy / maxima

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

Jacobi elliptic functions and complex arguments #2264

Open rtoy opened 3 months ago

rtoy commented 3 months ago

Imported from SourceForge on 2024-07-05 22:36:59 Created by vladutzanu on 2016-09-14 15:01:11 Original: https://sourceforge.net/p/maxima/bugs/3215


Hello

(This is a following from maxima-discuss: "Jacobi elliptic functions, maxima 5.38.1")

Given the following input:

kill(all)$ [wp, ws, Ap, As]:[1, 1.1, 1, 15]$
k  : wp / ws, numer$
k1 : (10^(Ap / 10) - 1) / (10^(As / 10) - 1), numer$
K  : elliptic_kc(k)$
K1 : elliptic_kc(k1)$
K_ : elliptic_kc(1-k)$
K1_: elliptic_kc(1-k1)$
N  : ceiling(K * K1_ / (K_ * K1)), numer;
lg : log(10), numer$
q0 : ((1 - sqrt(sqrt(1-k))) / (1 + sqrt(sqrt(1-k)))) / 2, numer$
q0sq : q0^4, numer$
q : q0 * (q0sq * (q0sq * (150 * q0sq + 15) + 2) + 1), numer$
As : 10 * log(0.0625 / q^N * (10^(Ap * 0.1) - 1) + 1) / lg, numer;
k1 : (10^(Ap * 0.1) - 1) / (10^(As * 0.1) - 1), numer$
q : (k1 / 16)^(1 / N), numer$
q_SQ : q^2, numer$
k : (q_SQ * (q_SQ * (q_SQ * (q_SQ * (q_SQ^2 + 2) + 2) + 1) + 2) + 1) /
    (q * (q * (q_SQ * (q * (4 * q_SQ * q + 8) + 4) + 4) + 4) + 1) * 4 * sqrt(q), numer$
K1 : elliptic_kc(k1), numer$
K  : elliptic_kc(k^2), numer$
R(x) := jacobi_cd(rectform(N * K1 / K * inverse_jacobi_cd(x / wp, k^2)), k1)$
H(x) := 1 / sqrt(1 + (10^(Ap / 10) - 1) * R(x)^2)$
R2(xi, x) := ((sqrt(1 - 1 / xi^2) + 1) * x^2 - 1) / ((sqrt(1 - 1 / xi^2) - 1) * x^2 + 1)$
p:1/k;
H2(x) := 1 / sqrt(1 + (10^(Ap / 10) - 1) * R2(R2(p, p), R2(p, x))^2)$
plot2d([H2(x), realpart(H(x)), imagpart(H(x))], [x, 0, 3], grid2d)$

R(x) can only be plotted with rectform() encompassing its argument, otherwise it ouputs these: BIGFLOAT: unable to convert 3.15105041004758 (2.257205326820854-6.664001874625061e-8%i) to a CL or BIGFLOAT number.

It should plot the frequency response of an elliptic/Cauer filter. The imagpart(H(x)) should be, ideally, zero.

rtoy commented 3 months ago

Imported from SourceForge on 2024-07-05 22:37:00 Created by robert_dodier on 2016-09-16 05:15:53 Original: https://sourceforge.net/p/maxima/bugs/3215/#9eaf


Thanks for the bug report. Can you please state specifically what is the expression which triggers the error. I can't guess what it is.

rtoy commented 3 months ago

Imported from SourceForge on 2024-07-05 22:37:03 Created by vladutzanu on 2016-09-17 05:47:20 Original: https://sourceforge.net/p/maxima/bugs/3215/#9eaf/bc90


R(x) is the culprit, 6th line from the bottom. If declared like this:

R(x) := jacobi_cd(N * K1 / K * inverse_jacobi_cd(x / wp, k^2), k1)$

It won't work, messages about BIGFLOAT and whatnot keep flooding the screen. The argument must be encompassed in rectform() to work, as it is now (a solution suggested in the mailing list). I think it's the inverse functions that behave abnormally.

rtoy commented 3 months ago

Imported from SourceForge on 2024-07-05 22:37:07 Created by robert_dodier on 2016-09-17 20:43:57 Original: https://sourceforge.net/p/maxima/bugs/3215/#9eaf/bc90/a753


OK, I get it now.

If I replace the definition of R(x) as you say, I get the same plot; I don't see any errors about BIGFLOAT.

I am working with:

Maxima version: "branch_5_38_base_223_gcf9cbb2"
Maxima build date: "2016-08-24 16:33:57"
Host type: "i686-pc-linux-gnu"
Lisp implementation type: "CLISP"
Lisp implementation version: "2.49 (2010-07-07) (built 3605610186) (memory 3681070445)"

Note that I made some bug fixes to the elliptical functions code on 2016-07-26 (commits 80cab10, 133408b, and 372582e). Maybe you can try your problem with current source code.

These bug fixes and others will appear in Maxima 5.39, which should be on the horizon. But I'm not making releases these days, so I don't know for sure.

rtoy commented 3 months ago

Imported from SourceForge on 2024-07-05 22:37:10 Created by robert_dodier on 2016-09-17 20:45:05 Original: https://sourceforge.net/p/maxima/bugs/3215/#8f6f


rtoy commented 3 months ago

Imported from SourceForge on 2024-07-05 22:37:14 Created by vladutzanu on 2016-09-20 20:00:46 Original: https://sourceforge.net/p/maxima/bugs/3215/#9eaf/bc90/a753/a58c


These bug fixes and others will appear in Maxima 5.39, which should be on the horizon. But I'm not making releases these days, so I don't know for sure.

If you say the changes were made, then that's the solution. I'll try to see how the changes cope with my current install (archlinux x64), if not, I'll wat, patiently, for the new release. The rectform() workaround will hold until then. Thank you for the fix.