CoCalc Public FilesTreball EDOS.sagewsOpen with one click!
Authors: Xavi Duran, Lara Reche Polo, Sara Roig Segura, Xènia Serra
Views : 78
#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())
x  43x2203x+253\displaystyle x \ {\mapsto}\ \frac{4}{3} \, x^{2} - \frac{20}{3} \, x + \frac{25}{3}
x  427x4+4027x3269x211027x+1175108\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}
x  42187x8+802187x75562187x6+13402187x5+18854374x469252187x3+201178748x2293358748x+1308575139968\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}
x  414348907x16+16014348907x159044782969x14+2492014348907x1312899914348907x12+1069904782969x11+20034728697814x10554492528697814x9+2177352551018336x819503625114791256x7237091877459165024x6+170265635153055008x590030365513673320192x4+56067402551836660096x321730010092448880128x2352263226157346640384x+2444639078975235092492288\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
[a(t)=2000e(1250t)+2000\displaystyle a\left(t\right) = -2000 \, e^{\left(-\frac{1}{250} \, t\right)} + 2000, b(t)=2000e(1250t)+1000e(1125t)+1000\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
[a1(t)=4687589591098454e(1250t)\displaystyle a_{1}\left(t\right) = \frac{468758959}{1098454} \, e^{\left(-\frac{1}{250} \, t\right)}, b1(t)=4687589591098454e(1250t)7737412012655592029662886230e(1125t)\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