lars b?ttcher <lars-boettcher at web.de> writes:
> Hi,
> I would like to solve an equation with maximas's solve-function. The problem
> is that it is an implicit function and the solve-function give no result. How can I
> solve this Problem.
>
>
> (%i1) eqn3: 1/w*sqrt(Rs(rho_c+0.2*Rs*h^2))*coth(sqrt(Rs/(rho_c+0.2*Rs*h^2))*d)-Rc;
>
> (%i2) Rs: 2.53;
>
> (%i3) Rc: 0.395;
>
> (%i4) h: 1.4e-6;
>
> (%i5) d: 49e-6;
>
> (%i6) w: 100e-6;
>
> (%i7) solve(eqn3,rho_c);
Differentiating eqn3, I convinced myself that there was only one root (I
think the function is increasing).
Then have another look at the equation. You'd be done if you could find
rho_c+0.2*Rs*h^2, which appears twice in the equation.
So let's substitute in a name for that:
(%i1) eqn3: 1/w*sqrt(Rs*(rho_c+Rs*h^2/5))*coth(sqrt(Rs/(rho_c+Rs*h^2/5))*d)-Rc;
2
h Rs Rs
sqrt(Rs (----- + rho_c)) coth(d sqrt(-------------))
5 2
h Rs
----- + rho_c
5
(%o1) ---------------------------------------------------- - Rc
w
(%i2) tau: sqrt(rho_c+Rs*h^2/5);
2
h Rs
(%o2) sqrt(----- + rho_c)
5
(%i3) ratsubst('tau^2, tau^2, eqn3);
d sqrt(Rs)
Rc w - sqrt(Rs) abs(tau) coth(----------)
abs(tau)
(%o3) - -----------------------------------------
w
This already looks rather more feasible. After a bit of rearrangement
(or just typing it in), you get that tau solves this equation if
d sqrt(Rs) Rc w
(%o7) abs(tau) coth(----------) = --------
abs(tau) sqrt(Rs)
I can't see a way to get an analytical solution here: exponentialising
doesn't seem to help because you end up with both exponentials from the
coth and the non-exponential terms from the abs(tau).
Now, since you're doing this with concrete numbers anyway, it doesn't
really matter if we numerically find a root. So, put in the values of
the constants (and lets get rid of the abs(tau) nonsense: I know that
tau is real and positive).
(%i20) subst([Rs = 2.53, d = 49e-6, w = 100e-6, Rc = 0.395], %o7);
7.793927123087563e-5
(%o20) abs(tau) coth(--------------------) = 2.4833437231746297e-5
abs(tau)
(%i21) subst(x,abs('tau),%);
7.793927123087563e-5
(%o21) coth(--------------------) x = 2.4833437231746297e-5
x
(here I realise that I shouldn't have actually bound a value to tau. Ho
hum, I'll just cheat and call it something else). I'll use the mnewton
package to find the root.
(%i27) load("mnewton")$
(%i28) mnewton([%o22],[x],[1]);
(%o28) [[x = 2.4742405969463985e-5]]
Then
(%i32) tau = rhs(first(first(%o28))), Rs=2.53, h=1.4e-6;
(%o32) sqrt(rho_c + 9.917599999999998e-13) = 2.4742405969463985e-5
and
(%i35) allroots(%o32^2);
(%o35) [rho_c = 6.111948931577672e-10]
Tada! (Note you don't use solve here, since that would replace our
inexact floating point values with rational approximations, which
doesn't make sense).
Exercise for reader: do it again with interval arithmetic and find
bounds for rho_c given bounds on the physical constants you're putting
in. I mean, how accurate is this answer?
Rupert
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 314 bytes
Desc: not available
Url : http://www.math.utexas.edu/pipermail/maxima/attachments/20090428/43677fe2/attachment.pgp