Commit 209627c1 authored by jben's avatar jben

good extraction of reals roots of a order 3 polyno

parent d46259e8
This diff is collapsed.
algebraic:true$
d_Ec:1/2*d_masse*(d_V1^2-d_V0^2)$
d_V1r:d_V1+d_vent$
d_V0r:d_V0+d_vent$
d_t:2*d_d/(d_V0+d_V1)$
d_Vrt:(t/d_t)*(d_V1r-d_V0r)+d_V0r$
d_Ea:1/2*d_rho_air*d_SCx*factor(integrate(d_Vrt^3,t,0,d_t))$
Poly:ratsimp((d_Ec+d_Ea+d_Ep+d_Er)/d_t-d_P,d_V1)$
a:ratcoef(Poly,d_V1^3)$
b:ratcoef(Poly,d_V1^2)$
c:ratcoef(Poly,d_V1)$
d:ev(Poly,d_V1=0)$
stringout("calc_model_a.out",ratsimp(a))$
stringout("calc_model_b.out",ratsimp(b))$
stringout("calc_model_c.out",ratsimp(c))$
stringout("calc_model_d.out",ratsimp(d))$
stringout("calc_model_Ea.out",ratsimp(d_Ea))$
stringout("calc_model_Ec.out",ratsimp(d_Ec))$
kill(a,b,c,d)$
p:(3*a*c-b**2)/(3*a^2)$
q:(2*b**3-9*a*b*c+27*a**2*d)/(27*a**3)$
ecart:-b/(3*a)$
stringout("calc_model_p.out",ratsimp(p))$
stringout("calc_model_q.out",ratsimp(q))$
kill(p,q)$
discri:27*q**2+4*p**3$
stringout("calc_model_discri.out",discri)$
Q:V**3+p*V+q$
forget(facts())$
assume(discri>0)$
sol:realpart(subst(solve(Q=0,V)[3],V))$
stringout("calc_model_sol_pos.out",ratsimp(sol)+ecart)$
forget(facts())$
assume(discri<0)$
sol:realpart(subst(solve(Q=0,V)[3],V))$
stringout("calc_model_sol_neg.out",ratsimp(sol)+ecart)$
......@@ -72,96 +72,66 @@ sub maxima2cc
return $out;
}
sub get_result
{
my $filename = $_[0];
my $varname = $_[1];
my $indent = $_[2];
open(my $g2,"<$filename");
my $inEa="";
while(<$g2>)
{
chomp;
$inEa="$inEa$_";
}
close($g2);
my $r=maxima2cc($inEa,$varname,$indent,"d_prov_${varname}_");
return $r;
}
# On écrit les equations
open(my $f,">calc_model.mc");
# Energie cinetique
print $f 'd_Ec:1/2*d_masse*(d_V1^2-d_V0^2)$';
print $f "\n";
# Energie depensé dans les frottement du vents (avec une progression
# affine de la vitesse entre V0 et V1
print $f 'd_V1r:d_V1+d_vent$';
print $f "\n";
print $f 'd_V0r:d_V0+d_vent$';
print $f "\n";
print $f 'd_Vr:x/d_d*(d_V1r-d_V0r)+d_V0r$';
print $f "\n";
print $f 'd_Ea:1/2*d_rho_air*d_SCx*factor(integrate(d_Vr^2,x,0,d_d))$';
print $f "\n";
# On calcule le temps en fonction des vitesses
print $f 't:2*d_d/(d_V0+d_V1)$';
print $f "\n";
# On écrit l'équation
print $f 'Equation:d_Ec+d_Ea+d_Ep+d_Er=t*d_P$';
print $f "\n";
# On la resoud et on récupere la racine
print $f 'r:subst(solve(Equation,d_V1)[3],d_V1)$';
print $f "\n";
# On ecrit les resultats
print $f 'stringout("calc_model_V1.out",r)$';
print $f 'stringout("calc_model_Ea.out",d_Ea)$';
print $f 'stringout("calc_model_Ec.out",d_Ec)$';
print $f "\n";
# On ferme le fichier
close($f);
# processing du calc_model.mc
system("maxima -b calc_model.mc");
# Lecture du resultat
open(my $g,"<calc_model_V1.out");
my $inV1="";
while(<$g>)
{
chomp;
$inV1="$inV1$_";
}
close($g);
open(my $g2,"<calc_model_Ea.out");
my $inEa="";
while(<$g2>)
{
chomp;
$inEa="$inEa$_";
}
close($g2);
open(my $g3,"<calc_model_Ec.out");
my $inEc="";
while(<$g3>)
{
chomp;
$inEc="$inEc$_";
}
close($g3);
# On ouvre la sortie
open(my $h,">calc_model.cc");
print $h "#include<cmath>\n\n";
print $h "void model::calc_model(double & d_V1, double & d_Ea, double & d_Ec, double d_V0, double d_vent, double d_d, double d_Ep, double d_Er) const\n";
print $h "{\n";
print $h " double a,b,c,d,p,q,discri;\n\n";
# on calcul a,b,c,d puis p,q puis discri
print $h get_result("calc_model_a.out","a"," ");
print $h get_result("calc_model_b.out","b"," ");
print $h get_result("calc_model_c.out","c"," ");
print $h get_result("calc_model_d.out","d"," ");
print $h get_result("calc_model_p.out","p"," ");
print $h get_result("calc_model_q.out","q"," ");
print $h get_result("calc_model_discri.out","discri"," ");
# puis si discri > 0 on prend une solution, sinon l'autre
print $h " if(discri>0)\n";
print $h " {\n";
print $h get_result("calc_model_sol_pos.out","d_V1"," ");
print $h " }\n";
print $h " else\n";
print $h " {\n";
print $h get_result("calc_model_sol_neg.out","d_V1"," ");
print $h " }\n";
print $h maxima2cc($inV1,"d_V1"," ","d_provV1_");
print $h " if(d_V1<0)\n";
print $h " d_V1=0;\n\n";
print $h " if(std::isnan(d_V1))\n";
print $h " d_V1=d_V0;\n\n";
# puis on calcule les energie
print $h maxima2cc($inEa,"d_Ea"," ","d_provEa_");
print $h maxima2cc($inEc,"d_Ec"," ","d_provEc_");
print $h get_result("calc_model_Ec.out","d_Ea"," ");
print $h get_result("calc_model_Ea.out","d_Ec"," ");
print $h "}\n";
close($h);
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment