223.html

Manual Maple 2.2.3.mws

Comandos Úteis e Simples

Seção 2.2: Resolvendo Equações

Session 2.2.3: dsolve

Fibonacci

> restart;

> fib := (-4-6*x-6*x^2)*Y(x)+(-1+x+x^2)^2*diff(Y(x),x,x);

fib := (-4-6*x-6*x^2)*Y(x)+(-1+x+x^2)^2*diff(Y(x),`$`(x,2))

> Order := 14;

> dsolve( {fib,Y(0)=1,D(Y)(0)=1}, Y(x), series );

> seq(combinat[fibonacci](n),n=1..14);

> dsolve( fib, Y(x) );

Y(x) = _C1/(-1+x+x^2)+_C2*(6*x^5+15*x^4-10*x^3-30*x^2+30*x)/(-6+6*x+6*x^2)

Combustion

> de := diff(y(x),x) = y(x)^2*(1-y(x));

> sn := dsolve({de,y(0)=10^(-k)},y(x));

sn := y(x) = 1/(LambertW(-1/exp(-1)*(-1+10^(-k))/(10^(-k))*exp(-(-1+10^(-k))/(10^(-k)))*exp(-x-1))+1)

> plot3d(rhs(sn),x=0..10^3,k=0..3,axes=BOXED);

[Maple Plot]

> trial := eval(sn,k=5);

trial := y(x) = 1/(LambertW(99999/exp(-1)*exp(99999)*exp(-x-1))+1)

> plot(rhs(trial),x=0.999e5..1.001e5);

[Maple Plot]

Lotka-Volterra

> restart;

> dsys := {diff(x(t),t)=x(t)*(1-y(t)),
diff(y(t),t)=s*y(t)*(x(t)-1) };

No texto, esta linha apresenta dsys[1] e dsys[2]; mas isso quebra uma das regras da "Limpeza da Área de trabalho", e o conjunto de dysys pode até estar em uma ordem diferente. Logo, nesta área de trabalho, o comando foi modificado para usar a substituição. Isso também fez o texto menor e mais claro.

> de := diff(x(y),y) = subs( x(t)=x(y), y(t)=y, subs(dsys,diff(x(t),t)/diff(y(t),t)) );

> _EnvAllSolutions := true;

> dsolve(de,x(y));

x(y) = exp((-s*LambertW(_NN1,-y^(-1/s)*exp(1/s*y+1/s*_C1-2*I/s*Pi*_Z2))+y+_C1-ln(y)-2*I*Pi*_Z2)/s)

> expand(%);

x(y) = 1/exp(LambertW(_NN1,-y^(-1/s)*exp(1/s*y+1/s*_C1-2*I/s*Pi*_Z2)))*exp(1/s*y)*exp(1/s*_C1)/(y^(1/s))/exp(1/s*Pi*_Z2*I)^2

> simplify(%);

x(y) = -LambertW(_NN1,-y^(-1/s)*exp(-(-y-_C1+2*I*Pi*_Z2)/s))

> subs(_Z2=0,%);

x(y) = -LambertW(_NN1,-y^(-1/s)*exp(-(-y-_C1)/s))

> x1 := subs(_NN1=0,%);

x1 := x(y) = -LambertW(0,-y^(-1/s)*exp(-(-y-_C1)/s))

> x2 := subs(_NN1=-1,%%);

x2 := x(y) = -LambertW(-1,-y^(-1/s)*exp(-(-y-_C1)/s))

> s := rand()/1.0e12;

> toplot := {[rhs(x1),y,y=0..10],[rhs(x2),y,y=0..10]};

$toplot := {[-LambertW(-1/y^3.668370381*exp(3.668370381*y+3.668370381*_C1)), y, y = 0 .. 10], [-LambertW(-1,-1/y^3.668370381*exp(3.668370381*y+3.668370381*_C1)), y, y = 0 .. 10]}$
$toplot := {[-LambertW(-1/y^3.668370381*exp(3.668370381*y+3.668370381*_C1)), y, y = 0 .. 10], [-LambertW(-1,-1/y^3.668370381*exp(3.668370381*y+3.668370381*_C1)), y, y = 0 .. 10]}$

> sols := `union`(seq(toplot,_C1={-7,-6,-5,-4,-3,-2,-3/2,-1-s})):

> plot(%,view=[0..20,0..10],colour=black);

[Maple Plot]

> eval([x1,x2],y=1);

> subs(_C1=-1-s,%);

> evalf(%);

> s := 's';

> solve(X1 = rhs(eval(x1,y=1)), _C1);

> subs(_Z5=0,%);

> subs(s=0,dsys);

> dsolve(%,{x(t),y(t)});

$[{y(t) = _C2}, {x(t) = _C1*exp(Int(1-y(t),t))}]$

> dsolve(dsys,{x(t),y(t)});

$[{x(t) = 0}, {y(t) = _C1*exp(-s*t)}], [{x(t) = RootOf(-Int(1/(_a*(LambertW(_NN2,exp(_a*s)/(_a^s)/exp(s*Pi*_Z5*I)^2*exp(_C1*s)*exp(-1))+1)),_a = `` .. _Z)+t+_C2)}, {y(t) = (-diff(x(t),t)+x(t))/x(t)}]$
$[{x(t) = 0}, {y(t) = _C1*exp(-s*t)}], [{x(t) = RootOf(-Int(1/(_a*(LambertW(_NN2,exp(_a*s)/(_a^s)/exp(s*Pi*_Z5*I)^2*exp(_C1*s)*exp(-1))+1)),_a = `` .. _Z)+t+_C2)}, {y(t) = (-diff(x(t),t)+x(t))/x(t)}]$

Airy

> restart;

> plots[setoptions](colour=BLACK);

> AiryDE := diff(y(z),z,z) = z*y(z);

AiryDE := diff(y(z),`$`(z,2)) = z*y(z)

> dsolve( AiryDE, y(z) ) ;

> HubbardWest := diff(x(t),t) = x(t)^2 - t;

> dsolve( HubbardWest, x(t) );

> dsolve( {HubbardWest, x(0)=alpha}, x(t) );

x(t) = -(-(3*GAMMA(2/3)^2*3^(2/3)+2*alpha*Pi*3^(5/6))/(-3*GAMMA(2/3)^2*3^(1/6)+2*alpha*Pi*3^(1/3))*AiryAi(1,t)+AiryBi(1,t))/(-(3*GAMMA(2/3)^2*3^(2/3)+2*alpha*Pi*3^(5/6))/(-3*GAMMA(2/3)^2*3^(1/6)+2*alph...

> X := rhs( % ):

> Nplot := 10;

>    P := plot( [seq(X,alpha=[seq(-5*i/Nplot, i=0..Nplot)])],
           t=0..8, y=-5..0,
           scaling=CONSTRAINED ): plots[display](P);

[Maple Plot]

> P1 := plot( eval(X,alpha=0), t=0..300, linestyle=4 ): plots[display](P1);

[Maple Plot]

> N := 3000;

> h := 300.0/N;

> badsol := dsolve({HubbardWest,x(0)=0},numeric,method=classical[rk4],stepsize=h,output=procedurelist);

> wrongpts := [seq( subs( badsol(k*h),[t,x(t)]), k=0..N)]:

> P2 := plot( wrongpts, style=POINT, colour=black, symbol=POINT ):

> plots[display]({P1,P2});

[Maple Plot]

> goodsol := dsolve( {HubbardWest,x(0)=0}, numeric, stiff=true, range=0..1.0e8 );

> plots[odeplot]( goodsol, style=POINT, symbol=POINT, colour=BLACK );

[Maple Plot]

> goodsol( 1.0e5 );

> eval( X, [alpha=0, t=1.0e5] );

-(-.1298219251e-9155730*sqrt(3)+.3876788305e9155733)/(.4105329705e-9155733*sqrt(3)+.1225948115e9155731)

> evalf( % );

> goodsol( 1.0e7 );

> -sqrt( 1.0e7 );

> eval( X, [alpha=0, t=1.0e7] );

> goodsol( 1.0e8 );

> -sqrt( 1.0e8 );

Caminhos diferentes para a mesma equação.

> restart;

>    Suzen := proc (m, t, yvec, ypvec)
    local i, j;
    option `[yvec[1] = x(t), yvec[2] = diff(x(t),t),
             yvec[3] = y(t), yvec[4] = diff(y(t),t)]`;
    for i to m/4 do
      j := 4*(i-1);
      ypvec[1+j] :=   yvec[2+j];
      ypvec[2+j] := -yvec[1+j]^2+yvec[3+j]^2;
      ypvec[3+j] :=   yvec[4+j];
      ypvec[4+j] := 2*yvec[1+j]*yvec[3+j];
    end do;
    ypvec
  end proc:

> N := 72;

> ic := map(evalf,vector( 4*N, [seq(op([cos(2*Pi*(i-1)/N),0,sin(2*Pi*(i-1)/N),0]),i=1..N)])):

> pvars := [seq(z||i(t),i=1..4*N)]:

> st := time(): sol := dsolve( numeric, procedure=Suzen, range=0..3, start=0, initial=ic, procvars=pvars ):
time_taken_de := time() - st;

Warning, cannot evaluate the solution further right of 2.9744778, probably a singularity

> st := time(): plots[odeplot]( sol, [seq([z||(1+4*(i-1))(t),z||(3+4*(i-1))(t)],i=1..N)], colour=black, view=[-3..3,-3..3], scaling=CONSTRAINED, labels=["",""], axes=BOXED );
time_taken_plotting := time() - st;

> time_taken_plotting/time_taken_de;

>