Plotar o Campo de Direções para um sistema de primeira ordem autônoma x' = f(x,y), y' = g(x,y) em um conjunto de N pontos (x1 , y1), ...,(xN,yN) dados, significa desenhar, para cada ponto (xJ,yJ), um vetor com origem neste ponto e com a direção do vetor V = (x',y') = ( f(x0,y0), g(x0,y0) ). O vetor V representa a "velocidade " com que a solução (x(t),y(t)) passa por este ponto.

O Maple V traça o Campo de Direções de um sistema de equações diferenciais através do comando dfieldplot da biblioteca DEtools

Campo := dfieldplot( [ D(x)(t) = f(x,y), D(y)(t) = g(x,y)], [x(t), y(t)], t=a..b, x=c..d, y=e..f , dirgrid[m,m], arrows=OPÇÃO, title = `titulo`, color= cor, e outras opções de plots);

Na opção "dirgrid" escolhemos uma grade de "N = m x n" pontos (o defaul é [20,20] ) e em "arrows=OPÇÃO" também podemos escolher THIN (seta fina), SLIM (seta delgada), THICK (seta espessa), LINE (segmento de reta) , NONE (nada), para o desenho das linhas tangentes.

Plotar o Retrato de Fase (ou Espaço de Fase) de um sistemea autônomo x' = f(x,y) , y' = g(x,y) em um conjunto em conjunto N pontos (x1 , y1), ...,(xN,yN) dados, significa plotar simultaneamente as N soluções ((x(t),y(t)) que passam por cada ponto (xJ,yJ) no instante , isto é, desenhar as soluções que satisfazem a condição inicial x(0)=xJ e y(0)=yJ

Com o comando phaseportrait da biblioteca DEtools o Maple V plota o Retrato de Fase da equação nos pontos e domínios especificados. Este comando plota automaticamente também o Campo de Direções. A sintaxe do comando é dada por

EspFase := phaseportrait( [ D(x)(t)=f(x,y), D(y)(t)=g(x,y)], [x(t), y(t)], t=a..b, [ [x(0)=a1, y(0)=b1] , ..., [x(0)=aM, y(0)=bM] ], x=c..d, y=e..f , dirgrid[m,m], arrows=OPÇÃO, title= `titulo`, color= cor, opções de dsolve/numeric, e outras opções de plot)

Vejamos alguns exemplos simples

> restart;

> with(DEtools):

> dfieldplot([D(x)(t)=x+y,D(y)(t)=y],[x(t),y(t)], t=-2..2, x=-5..5, y=-5..5, dirgrid=[10,10], arrows=LARGE);

> phaseportrait([D(x)(t)=x+y,D(y)(t)=y],[x(t),y(t)],t=-2..2, [[x(0)=1, y(0)=1]], linecolor=black);

Exemplo com algumas opções adicionais

> phaseportrait([D(x)(t)=x+y,\ D(y)(t)=y],[x(t),y(t)], t=-2..2,[[x(0)=1,y(0)=1]], x=-2..2,y=-2..2, dirgrid=[10,10], linecolor=black, color=x^2+y^2, method=classical, stepsize=0.1);

Exemplo da equação de Lotka-Volterra

> phaseportrait([diff(x(t),t)= x(t)*(1-y(t)), diff(y(t),t) =.3*y(t)*(x(t)-1)], [x(t),y(t)], t=0..20, [[x(0)=1, y(0)=1.5]], x=-1..2, y=-1..2, title=`Lotka-Volterra`, color = [.3*y(t)*(x(t)-1), x(t)*(1-y(t)), .1], linecolor=black);

>

Oscilações forçadas e amortecidas

Um corpo de massa M está preso a uma mola de constante elástica K e oscila em um meio viscoso tal que a força de atrito é proporcional à velocidade com constante L >0. Supomos que y(t) denota a coordenada do corpo no instante t. Vamos supor ainda que o corpo está sujeito a uma força externa periódica F(t) = A cos ( w t). A equação do movimento é dada por

M y''+ L y'+ K y =F(t)

Sabemos que a solução desta equação é dada por

y(t) = T(t) + R cos (w t - d)

na qual T(t), é solução da equação homogênea e .

O termo T(t) é chamado de transiente . A solução y(t) tende para R cos(w t - d) que é chamada de solução permanente . A amplitude é chamada de R e depende das constantes M, L, K, A e w .

Para efeito de simplificação vamos supor que M = K = 1, e, neste caso temos

É fácil ver que R(0) =1 e

Problema :

a) Estudar graficamente a variação de R(w) para os valores de L= 1.5, 1.2 , 1.0 e 0.25 e, em cada caso, determinar o valor aproximado da frequência w para a qual a amplitude R é máxima.

b) Obter a solução exata do problema e, em cada caso, comparar os gráficos da solução exata com a da solução aproximada

Solução

(a) Escrevendo o valor de R em função de w e L

> RLw := 1/( (1-w^2)^2+L^2*w^2 )^(1/2);

Vamos escrever os valores de L dados na forma de uma lista :

> Lis := 0.25, 1.0, 1.2, 1.5;

Vamos escrever uma seqüencia de funções R(w) para cada valor de L dados na lista acima

> Rw := seq( subs(L=Lis[i], RLw), i=1..4);

O gráfico de R(w) para os valores de L solicitados estão plotados abaixo. Clicando-se o mouse no ponto de máximo de cada gráfico obtemos os pontos e os valores de maximos aproximados

L = 0,25 L = 1,0 L = 1,2 L = 1,5

Se L = 0,25 então para w = 0,97 R(w) tem um valor máximo 4,0

Se L = 1,0 então para w = 0,72 R(w) tem um valor máximo de 1,2.

Se L =1,2 então para w = 0,52 R(w) tem um valor máximo de 1,05.

Se L = 1,5 então R(w ) é decrescente se , logo tem um máximo em w=0

As figuras acima podem ser obtidas com

> plot( Rw[1], w=0..2 );

> plot( Rw[2], w=0..2, y=0..2 );

> plot( Rw[3], w=0..2, y=0..2 );

> plot( Rw[4], w=0..2, y=0..2 );

(b) Vamos calcular a solução exata da equação amortecida forçada para L = 0,25.

Vamos trocar o coeficiente L=0.25, dado na forma de um float, para a forma fracionária L= 1/4 :

> OscAmort := diff(y(x),x$2)+(1/4)*diff(y(x),x)+ y(x) = cos(x);

> SolOsc := dsolve({OscAmort, y(0)=0, D(y)(0)=0}, y(x));

Vamos simplificar a resposta fornecida com combine(expressão,trig)

> SolOsc2 := rhs(combine(SolOsc, trig));

A primeira parcela de SolOsc2 (que apresenta a exp(-1/8 x) ) tende a zero quando x tende para infinito. A solução particular 4 sen(x) é a solução permanente.

Vamos comparar o gráfico da solução transiente Solosc2, com o o gráfico da solução permanente 4 sin(x) e depois plotar a diferença Sol2 - 4 sin(x)

> plot({SolOsc2-4*sin(x) }, x=0..30, y=-4..4):

> plot({4*sin(x) }, x=0..30):

> plot({SolOsc2 }, x=0..20, y=-4..4):

>

Sol Exata = SolOsc2 Sol. Peman.= 4*sin(x) Diferença = SolOsc2-4*sin(x)