solving strategy, getting float-numbers?
- Subject: solving strategy, getting float-numbers?
- From: Gloeckner, Robert
- Date: Tue, 10 Jan 2006 23:54:32 +0100
hello,
i try to use maxima for daily work, but somehow i do not know the right
solving strategies. maybe some of you can give me hints.
the following maxima-code works fairly, but even in the second solving
(with numbers) i am not able to get numbers for the vars, although they
only contain numbers themselves.
the first solving stops at solving for k7 - i left them in the code in
order to show what my intention was.
the second solving uses numerical data for the curves. i managed to get
all variables solved (somehow strange for k7 again), but i do not get
numbers, why?
thanks for infos,
robert
/* $Header: /home/cvs/lisp/test.max,v 1.10 2006/01/10 22:53:38 robert
Exp $
*
* (C) Robert Gloeckner, 2006
*
* This is a maxima-file about adaptive fitting pvT-data for
* 7(11)-coefficient-fit. the curve will be divided in two (three)
* parts: (1) a linear, an (2) exponential and (3) another linear part.
*
*
* the curve is described by 2 functions,
* for (1) and (2):
* v(k1, k2, k3, k4, k5, k6, k7) = vv(k1, k2, k3, k4) + vvv(k5, k6, k7)
*
* for (3)
* u(k8, k9, k10, k11)
*
*
* point 1->2 & 5->6 linear part 1 => vv(k1, k2, k3, k4)
* point 2->3 & 6->7 exponential part => vvv(k5, k6, k7)
* point 3->4 & 7->8 linear part 2 => u(k8, k9, k10, k11)
*
* so there are 8 points = 10 equations but 11 variables.
*
* i introduce a free parameter to have enough equations
* i hope i do not kill my solvability(?) (i fear so)
* the chosen free parameter describes the exponential fraction
* at point 2: L = vvv_2 / v_2
*/
kill(all);
/* linear part1 of the pvT-curve */
vv(p, T, k1, k2, k3, k4) := (k1 / (k3 + p)) + T * ( k2 / (k4 + p));
/* exponential part of the pvT-curve */
vvv(p, T, k5, k6, k7) := k7 * exp ( k5 * p + k6 * T);
/* 7-coefficient-model of the pvT-curve (linear part 1 + exponential
part) */
v(p, T, k1, k2, k3, k4, k5, k6, k7) := vv(p, T, k1, k2, k3, k4) + vvv(p,
T, k5, k6, k7);
/* linear part2 of the pvT-curve */
u(p, T, k8, k9, k10, k11) := vv(p, T, k8, k9, k10, k11);
/************** i n p u t - d a t a ***********************************/
v1 : v(p1, T1, k1, k2, k3, k4, k5, k6, k7) = v_1; /* curve-point 1 */
v2 : v(p1, T2, k1, k2, k3, k4, k5, k6, k7) = v_2; /* curve-point 2 */
v3 : v(p1, T3, k1, k2, k3, k4, k5, k6, k7) = v_3; /* curve-point 3 */
v3_: u(p1, T3, k8, k9, k10, k11) = v_3;
v4 : u(p1, T4, k8, k9, k10, k11) = v_4; /* curve-point 4 */
v5 : v(p2, T5, k1, k2, k3, k4, k5, k6, k7) = v_5; /* curve-point 5 */
v6 : v(p2, T6, k1, k2, k3, k4, k5, k6, k7) = v_6; /* curve-point 6 */
v7 : v(p2, T7, k1, k2, k3, k4, k5, k6, k7) = v_7; /* curve-point 7 */
v7_: u(p2, T7, k8, k9, k10, k11) = v_7;
v8 : u(p2, T8, k8, k9, k10, k11) = v_8; /* curve-point 8 */
v1a: vv(p1, T1, k1, k2, k3, k4) = v_1; /* approximation vvv(T1) -> 0
*/
v5a: vv(p2, T5, k1, k2, k3, k4) = v_5; /* approximation vvv(T5) -> 0
*/
v2a: vvv(p1, T2, k5, k6, k7) = L * v_2; /* exponential fraction in
point 2 as _FREE_PARAMETER_ */
v2b: vv(p1, T2, k1, k2, k3, k4) = ( 1 - L) * v_2; /* resulting linear
fraction in point 2 */
/************ s y m b o l i c - s o l v i n g ********************/
kill( k1, k2, k3, k4, k5, k6, k7, k8, k9, k10, k11);
solveradcan : true;
globalsolve : true;
/* */
k7 : part( part( solve( ev( v2a), k7), 1), 2);
/* k5 v2a */
k1 : part( part( solve( ev( v2b), k1), 1), 2);
/* k1 k5 v2a v2b */
k2 : part( part( solve( ev( v5a), k2), 1), 2);
k1 : ev( k1)$ k7 : ev( k7)$ k2 : ev( k2)$
/* k1 k2 k7 v2a v2b v5a */
k3 : part( part( solve( ev( v1a), k3), 1), 2);
k1 : ev( k1)$ k7 : ev( k7)$ k2 : ev( k2)$ k3 : ev( k3)$
/* k1 k2 k3 k7 v1a v2a v2b v5a */
k6 : part( part( solve( ev( v3), k6), 1), 2);
k1 : ev( k1)$ k7 : ev( k7)$ k2 : ev( k2)$ k3 : ev( k3)$ k6 : ev(
k6)$
/* k1 k2 k3 k4 k7 v1a v2a v2b v3 v5a */
k4 : part( part( solve( ev( v6), k4), 1), 2);
k1 : ev( k1)$ k7 : ev( k7)$ k2 : ev( k2)$ k3 : ev( k3)$ k6 : ev( k6)$ k4
: ev( k4)$
/* k1 k2 k3 k4 k7 k6 v1a v2a v2b v3 v5a v6 */
t: ratsimp( ev( v7));
k7: rhs(solve( t, k7)[1]);
k1:ev(k1)$ k7:ev(k7)$ k2:ev(k2)$ k3:ev(k3)$ k6:ev(k6)$ k4:ev(k4)$
k5:ev(k5)$
for f in ['k1, 'k2, 'k3, 'k4, 'k5, 'k6, 'k7] do print( f = ev( f,
numer));
/**********e x a m p l e - d a t a **************************/
T5 : T1;
L : 1e-2;
T1 : 100; T2 : 200; T3 : 220; T4 : 300;
T5 : 100; T6 : 210; T7 : 240; T8 : 330;
p1 : 1; p2 : 1000;
v_1 : 20; v_2 : 40; v_3 : 55; v_4 : 80;
v_5 : 10; v_6 : 20; v_7 : 40; v_8 : 60;
/************ ( n u m e r i c a l) - s o l v i n g ************/
kill( k1, k2, k3, k4, k5, k6, k7, k8, k9, k10, k11);
solveradcan : true;
globalsolve : true;
/* */
k7 : part( part( solve( ev( v2a), k7), 1), 2);
/* k5 v2a */
k1 : part( part( solve( ev( v2b), k1), 1), 2);
/* k1 k5 v2a v2b */
k2 : part( part( solve( ev( v5a), k2), 1), 2);
k1 : ev( k1)$ k7 : ev( k7)$ k2 : ev( k2)$
/* k1 k2 k7 v2a v2b v5a */
k3 : part( part( solve( ev( v1a), k3), 1), 2);
k1 : ev( k1)$ k7 : ev( k7)$ k2 : ev( k2)$ k3 : ev( k3)$
/* k1 k2 k3 k7 v1a v2a v2b v5a */
k6 : part( part( solve( ev( v3), k6), 1), 2);
k1 : ev( k1)$ k7 : ev( k7)$ k2 : ev( k2)$ k3 : ev( k3)$ k6 : ev(
k6)$
/* k1 k2 k3 k4 k7 v1a v2a v2b v3 v5a */
k4 : part( part( solve( ev( v6), k4), 1), 2);
k1 : ev( k1)$ k7 : ev( k7)$ k2 : ev( k2)$ k3 : ev( k3)$ k6 : ev( k6)$ k4
: ev( k4)$
/* k1 k2 k3 k4 k7 k6 v1a v2a v2b v3 v5a v6 */
t: ratsimp( ev( v7));
t1: subst( exp(999*k5), kk5, solve( subst( kk5, exp(999*k5), t), kk5)) ;
k5: rhs(solve( t1^(1/999), k5)[1]);
k1:ev(k1)$ k7:ev(k7)$ k2:ev(k2)$ k3:ev(k3)$ k6:ev(k6)$ k4:ev(k4)$
k5:ev(k5)$
for f in ['k1, 'k2, 'k3, 'k4, 'k5, 'k6, 'k7] do print( f = ev( f,
numer));
--
Dipl.-Phys. Robert Gloeckner
Research Assistant
Deutsches Kunststoff-Institut DKI
German Institute for Polymers
Department of Technology
Schlossgartenstr. 6
D-64289 Darmstadt, Germany
phone +49(0)6151 - 16 6516
fax +49(0)6151 - 29 2855
http://www.DKI-online.de
GnuPG-Key: 1024D/9A0A2D72