On 5/28/10 3:31 AM, Je Suis wrote:
> Hello
>
> I did a search for "elliptic", "solve" and a few other keywords, but couldn't find anything related to my problem so, in case this isn't a trivial task that slips my eyes, I'll ask for your help:
>
> I have a Cauer filter where I have to calculate ws. This implies the use of elliptic_kc(k), where k has two forms: k2=1/ws^2 and k2`=sqrt(1-1/ws^2), so there are two integrals: K2 and K2`. These two make an equation something like this: C=K2`/K2, where C is known. From this I have to extract ws and I can't. For the record, if my mathematical skills would have been better, I wouldn't have been needed extra software... Here's what I did:
>
> --
> (some values to help)
> N:4$ wp:1$ Ap:3$ As:40$
> ep:sqrt(10^(Ap/10)-1)$
> es:1/sqrt(10^(As/10)-1)$
> t1:ep*es$
> t2:wp/ws$
> t1_:sqrt(1-t1^2)$
> t2_:sqrt(1-t2^2)$
>
> K1:elliptic_kc(t1^2)$ float(%); //reason...
> K1_:elliptic_kc(t1_^2)$ float(%); //...later
> K2:integrate(1/sqrt(1-t2^2),arg,0,%pi/2)$ float(%);
> K2_:integrate(1/sqrt(1-t2_^2),arg,0,%pi/2)$ float(%);
t2 is independent of arg, so I don't really follow what you're trying to
do here.
> --
>
> First, I tried the direct approach, with solve(k=elliptic_kc(t2_)/elliptic_kc(t2),ws), but I got this in return: elliptic_kc((ws^2-1)/ws^2)=elliptic_kc(1/ws^2)*k -- not helpful.
Yeah, maxima doesn't know anything about solving this. But note that
maxima knows make_ellipti_f(inverse_jacobi_sn(x, m)) is
elliptic_f(asin(x),m) and elliptic_f(z,%pi/2) = elliptic_kc(z). (Maxima
doesn't seem to know this latter equation.)
So you could express everything in terms of inverse_jacobi_sn and then
take the jacobi_sn of both sides. Don't know if that helps or not.
>
> Then I tried writing with integral form:
> solve(K=integrate((1/sqrt(1-(1-1/ws^2)*sin(phi)^2)),phi,0,%pi/2)/integrate((1/sqrt(1-1/ws^2*sin(phi)^2)),phi,0,%pi/2),ws)
> but I got this, in return:
> integrate(1/sqrt((1-sin(phi)^2)*ws^2+sin(phi)^2),phi,0,%pi/2)=integrate(1/(sqrt(ws-sin(phi))*sqrt(ws+sin(phi))),phi,0,%pi/2)*K - still not good.
>
> I then tried to ease up things by replacing k=N*K1/K1_=1.04829 -- no luck. Then I tried to "cheat" by trying to calculate the integral with the value ws should have, 1.3461 (where I had a bad surprise; I'll get to this later). I had to replace sin(phi) with 1 and use elliptic_kc(t1_) for K1_ -- see later on why -- and I got these values:
>
> K2:integrate(1/sqrt(1-t2^2),phi,0,%pi/2)$ float(%);
> K2_:integrate(1/sqrt(1-t2_^2),phi,0,%pi/2)$ float(%);
> 1.5708/sqrt(1.0-1.0/ws^2)
> 1.5708*abs(ws)
>
> With these I tried this:
>
> allroots(solve(float(K2/K2_=k),ws)) =>
> allroots: expected a polynomial; found [abs(ws)=(10479*abs(ws))/(10985*sqrt(ws^2-1))]
>
> Good enough. I simplified by hand and got: ws=sqrt(1/k^2+1)=1.38203. Close, but not quite, so I "blame" the elliminating of sin(phi). At this point, I decided to check the results with the known value of ws. Now, when I did the comparisons, I eliminated "sin(phi)^2", like this:
>
> integrate(1/sqrt(1-1.3461^2),phi,0,%pi/2)
>
> since it wouldn't calculate with a "sin(phi)" and replacing phi with %pi/2 would have been one. The results were the expected ones _except_ one: K1_=157.445 (result of integration) instead of 5.99391 (value given by elliptic_kc(t1_)). I then thought to install Scilab and make a comparison. I use Ubuntu 10.04 x86_64, installed, ran this:
>
> K1_=integrate('1/sqrt(1-t1_^2*sin(phi)^2)','phi',0,%pi/2)
> K1_=5.9939133
>
> Then I thought I would verify the error and did this:
>
> K1_=integrate('1/sqrt(1-t1_^2)','phi',0,%pi/2)
> K1_=157.44518
>
> which is a confirmation to me.
>
> If you managed to read this far, please let me know how can I calculate the integral(s) _with_ sin(phi) so I can use solve() to find out ws. It seems to work if I have numbers and Scilab, to my knowledge, isn't a symbolic solver. All I want is a formula to compute ws, that's all... Also, all these will be implemented in a spice program to calculate "on the spot" the poles and zeros, program which doesn't know what elliptic_kc() and such mean, it only has a handful of mathematical formulae, so approximations will be used instead. For the record, inverse Chebyshev is implemented, that means sin, acosh, etc.
Unfortunately, maxima doesn't know that
integrate(1/sqrt(1-m*sin(phi)^2), phi) is an elliptic integral, but
maxima does know some properties of elliptic_f and inverse_jacobi_sn.
Ray