Treball EDOS.sagews

#EXERCICI 1
##APARTAT a)
t=var("t")
s=var("s")
a=var("a")
z=var("z")
x=function("x",s)
y=function("y",s)
##Per poder trobar solucions generals del sistema d'equacions diferencials, calculem el polinomi característic de la matriu del sistema:
A=matrix([[a,-1],[3*a+13,1]])
pol=charpoly(A,z)
solve(pol(z),z)
##Els casos que tindrem dependran de com siguin els valors propis del polinomi característic.

###Cas (a+3)*(a-17)>0
eqx=diff(x(s),s)==(a*x(s)-y(s))
eqy=diff(y(s),s)==((3*a+13)*x(s)+y(s))

###El cas en què (a-17) i (a+3) siguin positius no dóna cap resultat. La suposició és inconsistent. Per això, agafem un cas particular, per exemple, que a=18.
eqx1=diff(x(s),s)==(18*x(s)-y(s))
eqy1=diff(y(s),s)==((3*18+13)*x(s)+y(s))
sols1=desolve_system([eqx1,eqy1],[x,y],ivar=s)
sols1[0].rhs().subs(s=((4*t^3)/3)+t);
sols1[1].rhs().subs(s=((4*t^3)/3)+t);

##El cas en què (a-17) i (a+3) siguin negatius no dóna cap resultat. La suposició és inconsistent. Per això, agafem un cas particular, per exemple, que a=-5.
eqx2=diff(x(s),s)==(-5*x(s)-y(s))
eqy2=diff(y(s),s)==((3*(-5)+13)*x(s)+y(s))
sols2=desolve_system([eqx2,eqy2],[x,y],ivar=s)
sols2[0].rhs().subs(s=((4*t^3)/3)+t);
sols2[1].rhs().subs(s=((4*t^3)/3)+t);

###Cas (a+3)*(a-17)<0
assume((a+3)>0)
assume((a-17)<0)
sols3=desolve_system([eqx,eqy],[x,y],ivar=s)
sols3[0].rhs().subs(s=((4*t^3)/3)+t);
sols3[1].rhs().subs(s=((4*t^3)/3)+t);

###Cas a=17
eqx4=diff(x(s),s)==(17*x(s)-y(s))
eqy4=diff(y(s),s)==((3*17+13)*x(s)+y(s))
sols4=desolve_system([eqx4,eqy4],[x,y],ivar=s)
sols4[0].rhs().subs(s=((4*t^3)/3)+t);
sols4[1].rhs().subs(s=((4*t^3)/3)+t);

[z == 1/2*a - 1/2*sqrt(a^2 - 14*a - 51) + 1/2, z == 1/2*a + 1/2*sqrt(a^2 - 14*a - 51) + 1/2]

/cocalc/lib/python2.7/site-packages/smc_sagews/sage_server.py:1013: DeprecationWarning: Substitution using function-call syntax and unnamed arguments is deprecated and will be removed from a future release of Sage; you can use named arguments instead, like EXPR(x=..., y=...) See http://trac.sagemath.org/5930 for details. exec compile(block+'\n', '', 'single') in namespace, locals

1/21*(sqrt(21)*(17*x(0) - 2*y(0))*sinh(1/6*sqrt(21)*(4*t^3 + 3*t)) + 21*cosh(1/6*sqrt(21)*(4*t^3 + 3*t))*x(0))*e^(38/3*t^3 + 19/2*t) 1/21*(sqrt(21)*(134*x(0) - 17*y(0))*sinh(1/6*sqrt(21)*(4*t^3 + 3*t)) + 21*cosh(1/6*sqrt(21)*(4*t^3 + 3*t))*y(0))*e^(38/3*t^3 + 19/2*t) -1/11*(sqrt(11)*(3*x(0) + y(0))*sinh(1/3*sqrt(11)*(4*t^3 + 3*t)) - 11*cosh(1/3*sqrt(11)*(4*t^3 + 3*t))*x(0))*e^(-8/3*t^3 - 2*t) -1/11*(sqrt(11)*(2*x(0) - 3*y(0))*sinh(1/3*sqrt(11)*(4*t^3 + 3*t)) - 11*cosh(1/3*sqrt(11)*(4*t^3 + 3*t))*y(0))*e^(-8/3*t^3 - 2*t) (cos(1/6*(4*t^3 + 3*t)*sqrt(-a^2 + 14*a + 51))*x(0) + ((a + 1)*x(0) - 2*x(0) - 2*y(0))*sin(1/6*(4*t^3 + 3*t)*sqrt(-a^2 + 14*a + 51))/sqrt(-a^2 + 14*a + 51))*e^(1/6*(4*t^3 + 3*t)*(a + 1)) (cos(1/6*(4*t^3 + 3*t)*sqrt(-a^2 + 14*a + 51))*y(0) + (2*a*(3*x(0) - y(0)) + (a + 1)*y(0) + 26*x(0))*sin(1/6*(4*t^3 + 3*t)*sqrt(-a^2 + 14*a + 51))/sqrt(-a^2 + 14*a + 51))*e^(1/6*(4*t^3 + 3*t)*(a + 1)) 8/3*(4*t^3 + 3*t)*e^(12*t^3 + 9*t)*x(0) - 1/3*(4*t^3 + 3*t)*e^(12*t^3 + 9*t)*y(0) + e^(12*t^3 + 9*t)*x(0) 64/3*(4*t^3 + 3*t)*e^(12*t^3 + 9*t)*x(0) - 8/3*(4*t^3 + 3*t)*e^(12*t^3 + 9*t)*y(0) + e^(12*t^3 + 9*t)*y(0)

##Ara passem a fer els casos particulars.
t=var("t")
s=var("s")
x=function("x",s)
y=function("y",s)
####Cas a=-1, amb condicions inicials t0=0, x0=1, y0=0
eqx5=diff(x(s),s)==((-1)*x(s)-y(s))
eqy5=diff(y(s),s)==((3*(-1)+13)*x(s)+y(s))
sols5=desolve_system([eqx5,eqy5],[x,y],ics=[0,1,0],ivar=s)
Sx=sols5[0].rhs().subs(s=((4*t^3)/3)+t);
Sy5=sols5[1].rhs().subs(s=((4*t^3)/3)+t);
p5=parametric_plot([Sx5,Sy5], (t,-1,1))
####Cas a=-1, amb condicions inicials t0=0, x0=2, y0=0
sols5bis=desolve_system([eqx5,eqy5],[x,y],ics=[0,2,0],ivar=s)
Sx5bis=sols5bis[0].rhs().subs(s=((4*t^3)/3)+t);
Sy5bis=sols5bis[1].rhs().subs(s=((4*t^3)/3)+t);
p5bis=parametric_plot([Sx5bis,Sy5bis], (t,-1.5,1.5))
show(p5+p5bis)

####Cas a=-2, amb condicions inicials t0=0, x0=1, y0=0
eqx6=diff(x(s),s)==((-2)*x(s)-y(s))
eqy6=diff(y(s),s)==((3*(-2)+13)*x(s)+y(s))
sols6=desolve_system([eqx6,eqy6],[x,y],ics=[0,1,0],ivar=s)
Sx6=sols6[0].rhs().subs(s=((4*t^3)/3)+t);
Sy6=sols6[1].rhs().subs(s=((4*t^3)/3)+t);
p6=parametric_plot([Sx6,Sy6], (t,-1.5,1.5))
####Cas a=-2, amb condicions inicials t0=0, x0=0.5, y0=0.5
sols6bis=desolve_system([eqx6,eqy6],[x,y],ics=[0,0.5,0.5],ivar=s)
Sx6bis=sols6bis[0].rhs().subs(s=((4*t^3)/3)+t);
Sy6bis=sols6bis[1].rhs().subs(s=((4*t^3)/3)+t);
p6bis=parametric_plot([Sx6bis,Sy6bis], (t,-1.5,1.5))
show(p6+p6bis)

####Cas a=-3, amb condicions inicials t0=0, x0=1, y0=0
eqx7=diff(x(s),s)==((-3)*x(s)-y(s))
eqy7=diff(y(s),s)==((3*(-3)+13)*x(s)+y(s))
sols7=desolve_system([eqx7,eqy7],[x,y],ics=[0,1,0],ivar=s)
Sx7=sols7[0].rhs().subs(s=((4*t^3)/3)+t);
Sy7=sols7[1].rhs().subs(s=((4*t^3)/3)+t);
p7=parametric_plot([Sx7,Sy7], (t,-0.75,0.75))
####Cas a=-3, amb condicions inicials t0=0, x0=-1, y0=0
sols7bis=desolve_system([eqx7,eqy7],[x,y],ics=[0,-1,0],ivar=s)
Sx7bis=sols7bis[0].rhs().subs(s=((4*t^3)/3)+t);
Sy7bis=sols7bis[1].rhs().subs(s=((4*t^3)/3)+t);
p7bis=parametric_plot([Sx7bis,Sy7bis], (t,-0.75,0.75))
show(p7+p7bis)

####Cas a=-4, amb condicions inicials t0=0, x0=1, y0=0
eqx8=diff(x(s),s)==((-4)*x(s)-y(s))
eqy8=diff(y(s),s)==((3*(-4)+13)*x(s)+y(s))
sols8=desolve_system([eqx8,eqy8],[x,y],ics=[0,1,0],ivar=s)
Sx8=sols8[0].rhs().subs(s=((4*t^3)/3)+t);
Sy8=sols8[1].rhs().subs(s=((4*t^3)/3)+t);
p8=parametric_plot([Sx8,Sy8], (t,-0.25,1.25))
####Cas a=-4, amb condicions inicials t0=0, x0=-1, y0=0
sols8bis=desolve_system([eqx8,eqy8],[x,y],ics=[0,-1,0],ivar=s)
Sx8bis=sols8bis[0].rhs().subs(s=((4*t^3)/3)+t);
Sy8bis=sols8bis[1].rhs().subs(s=((4*t^3)/3)+t);
p8bis=parametric_plot([Sx8bis,Sy8bis], (t,-0.25,1.25))
show(p8+p8bis)

Error in lines 3-3 Traceback (most recent call last): File "/cocalc/lib/python2.7/site-packages/smc_sagews/sage_server.py", line 1234, in execute flags=compile_flags), namespace, locals) File "", line 1, in <module> File "sage/calculus/var.pyx", line 133, in sage.calculus.var.function (build/cythonized/sage/calculus/var.c:1608) def function(s, **kwds): TypeError: function() takes exactly 1 positional argument (2 given)

##APARTAT b)
###Les solucions del sistema seran periòdiques si el discriminant que apareix en les solucions del polinomi característic és negatiu, i a més a més, si la part real dels valors propis és zero. El discriminant és positiu a (-inf,-3) i a (17,+inf), és zero a a=-3 i a=17, i és negatiu a (-3,17). Per tant, els possibles valors de "a" que busquem són els que es troben a l'interval (-3,17). La part real dels valors propis és (a+1)/2, per tant, això serà zero si i només si a=-1.

B=matrix([[-1,-1],[10,1]])
B
###La traça és nul·la, i el determinant és positiu. Mirant el pla traça-determinant, això correspon a tenir centres, és a dir, solucions periòdiques. La gràfica d'aquestes solucions es pot veure en l'apartat a).

[-1 -1] [10 1]

#EXERCICI 2
##Creem una funció que ens calculi la successió de polinomis a partir del polinomi inicial, la funció i l'iterat N.
def polinomi(polinicial, funcio, N):
    if N==0:
        return(polinicial)
    if N>=1:
        p0=polinicial;
        var("i,t,s");
        for i in range(1,N+1):
            p1=p0+(1/2)*(funcio-p0^2)
            p0=p1
        return(p0)

##Calculem els 7 primers polinomis. Per això, declarem 8 variables que ens guardin els polinomis, de manera que després sigui més senzill
##representar-los en un plot.
x=var("x")
f=function("f",x)
f(x)=sqrt((10-4*x)^2)
p0=function("p0",x)
p1=function("p1",x)
p2=function("p2",x)
p3=function("p3",x)
p4=function("p4",x)
p5=function("p5",x)
p6=function("p6",x)
p7=function("p7",x)
p0(x)=6*polinomi(0,((10-4*x)^2)/36,0)
p1(x)=6*polinomi(0,((10-4*x)^2)/36,1)
p2(x)=6*polinomi(0,((10-4*x)^2)/36,2)
p3(x)=6*polinomi(0,((10-4*x)^2)/36,3)
p4(x)=6*polinomi(0,((10-4*x)^2)/36,4)
p5(x)=6*polinomi(0,((10-4*x)^2)/36,5)
p6(x)=6*polinomi(0,((10-4*x)^2)/36,6)
p7(x)=6*polinomi(0,((10-4*x)^2)/36,7)
plot([f(x),p0(x),p1(x),p2(x),p3(x),p4(x), p5(x),p6(x),p7(x)], (x,0,5), ymin=-1, ymax=10, color=["red","greenyellow","forestgreen", "seagreen","lightseagreen","lightskyblue", "thistle", "mediumorchid"], legend_label=["f(x)","p0","p1","p2","p3","p4","p5","p6","p7"])
show(p1.expand())
show(p2.expand())
show(p3.expand())
show(p4.expand())

$\displaystyle x \ {\mapsto}\ \frac{4}{3} \, x^{2} - \frac{20}{3} \, x + \frac{25}{3}$

$\displaystyle x \ {\mapsto}\ -\frac{4}{27} \, x^{4} + \frac{40}{27} \, x^{3} - \frac{26}{9} \, x^{2} - \frac{110}{27} \, x + \frac{1175}{108}$

$\displaystyle x \ {\mapsto}\ -\frac{4}{2187} \, x^{8} + \frac{80}{2187} \, x^{7} - \frac{556}{2187} \, x^{6} + \frac{1340}{2187} \, x^{5} + \frac{1885}{4374} \, x^{4} - \frac{6925}{2187} \, x^{3} + \frac{20117}{8748} \, x^{2} - \frac{29335}{8748} \, x + \frac{1308575}{139968}$

$\displaystyle x \ {\mapsto}\ -\frac{4}{14348907} \, x^{16} + \frac{160}{14348907} \, x^{15} - \frac{904}{4782969} \, x^{14} + \frac{24920}{14348907} \, x^{13} - \frac{128999}{14348907} \, x^{12} + \frac{106990}{4782969} \, x^{11} + \frac{200347}{28697814} \, x^{10} - \frac{5544925}{28697814} \, x^{9} + \frac{21773525}{51018336} \, x^{8} - \frac{19503625}{114791256} \, x^{7} - \frac{237091877}{459165024} \, x^{6} + \frac{170265635}{153055008} \, x^{5} - \frac{9003036551}{3673320192} \, x^{4} + \frac{5606740255}{1836660096} \, x^{3} - \frac{2173001009}{2448880128} \, x^{2} - \frac{35226322615}{7346640384} \, x + \frac{2444639078975}{235092492288}$

#EXERCICI 3
##APARTAT a)
##Escrivim el sistema d'equacions diferencials que descriu el problema i el resolem.
var('t');
a=function("a",t);
b=function("b",t);
eqa=diff(a,t)==8-((4*a(t))/1000);
eqb=diff(b,t)==((4*a(t))/1000)-((8*b(t))/1000);
sols3=desolve_system([eqa,eqb],[a,b],ics=[0,0,0]);
show(sols3)

##Substituïm la variable t per 60 minuts = 1 hora.
sols3[0].rhs().subs(t=60);_.n()
sols3[1].rhs().subs(t=60);_.n()

#En el dipòpsit A hi ha una quantitat de 426.744277866893 g i en el dipòsit B hi ha una quantitat de 45.5276696730339 g. 

t

[$\displaystyle a\left(t\right) = -2000 \, e^{\left(-\frac{1}{250} \, t\right)} + 2000$, $\displaystyle b\left(t\right) = -2000 \, e^{\left(-\frac{1}{250} \, t\right)} + 1000 \, e^{\left(-\frac{1}{125} \, t\right)} + 1000$]

-2000*e^(-6/25) + 2000 426.744277866893 -2000*e^(-6/25) + 1000*e^(-12/25) + 1000 45.5276696730339

##APARTAT b)
##Escribim el nou sistema utilitzant la solució de l'apartat anterior com a condició inicial i el resolem:
var('t');
a1=function("a1",t);
b1=function("b1",t);
eqa1=diff(a1,t)==-((4*a1(t))/1000);
eqb1=diff(b1,t)==((4*a1(t))/1000)-((8*b1(t))/1000);
sols4=desolve_system([eqa1,eqb1],[a1,b1],ics=[0,426.744277866893,45.5276696730339]);
show(sols4)

t

[$\displaystyle a_{1}\left(t\right) = \frac{468758959}{1098454} \, e^{\left(-\frac{1}{250} \, t\right)}$, $\displaystyle b_{1}\left(t\right) = \frac{468758959}{1098454} \, e^{\left(-\frac{1}{250} \, t\right)} - \frac{773741201265559}{2029662886230} \, e^{\left(-\frac{1}{125} \, t\right)}$]

x=var('x')
solA=sols4[0].rhs();
solA
solB=sols4[1].rhs().subs(x=e^(-t/250));
solB

468758959/1098454*e^(-1/250*t) 468758959/1098454*e^(-1/250*t) - 773741201265559/2029662886230*e^(-1/125*t)

##Resolem les equacions solA=1 i solB=1 per tal de trobar el temps desitjat:
###solA=1
solve(solA==1,t)
solve(468758959/1098454*exp(-1/250*t) - 773741201265559/2029662886230*exp(-1/125*t)==1,t)

[t == 250*log(468758959/1098454)] [e^(1/125*t) == -773741201265559/1847745*e^(1/250*t)/(1098454*e^(1/250*t) - 468758959)]

##solB=1
solucio=solve(468758959/1098454*x - 773741201265559/2029662886230*x^2==1,x)
arrel1=solucio[0].rhs()
arrel1
solve(e^(-t/250)==arrel1,t)
250*log(-1547482402531118/3/(sqrt(82658769969849193301102896305) - 288715674232485)).n()

-3/1547482402531118*sqrt(82658769969849193301102896305) + 866147022697455/1547482402531118 [t == 250*log(-1547482402531118/3/(sqrt(82658769969849193301102896305) - 288715674232485))] 1513.52125636390

arrel2=solucio[1].rhs()
arrel2
solve(e^(-t/250)==arrel2,t)
n(250*log(1547482402531118/3/(sqrt(82658769969849193301102896305) + 288715674232485)))

3/1547482402531118*sqrt(82658769969849193301102896305) + 866147022697455/1547482402531118 [t == 250*log(1547482402531118/3/(sqrt(82658769969849193301102896305) + 288715674232485))] -27.6793216087572