MS 211 – Calculo Numérico
Projeto
Data de entrega: 16/06/2015
Observações para apresentação do projeto
• A nota do projeto comporá a média final.
• O projeto pode ser desenvolvido em grupo, com no máximo dois alunos. Não será permitido projetos com
três ou mais alunos.
• O projeto deverá ser entregue no formato impresso, e não será aceito após a data de entrega 16/06.
• As rotinas computacionais utilizadas também deverão ser anexadas ao projeto no formato impresso, independentemente da linguagem de programação utilizada.
• A nota do projeto será composta da avaliação dos resultados apresentados (teóricos e computacionais).
Soluções “soltas” ou “desconexas” não serão consideradas.
• Todos os métodos numéricos utilizados e que estejam cobertos pela ementa da disciplina deverão ser implementados.
Sugestões para elaboração do projeto
• Justifique todas as etapas do desenvolvimento do projeto.
• Faça um resumo dos métodos numéricos utilizados.
• Utilize recursos gráficos para auxiliar o desenvolvimento do projeto.
• Apresente os algoritmos dos métodos numéricos utilizados.
• Inclua as listagens dos programas.
• Apresente claramente os resultados numéricos obtidos fazendo uma análise objetiva.
• Apresente conclusão e bibliografia.
1
Atividade 01
Um pêndulo simples é um sistema composto de uma massa pontual acoplada em um fio de massa desprezível e de
comprimento L. O fio está preso a uma extremidade permitindo um movimento oscilatório da massa. Quando a
gravidade é a única força agindo sobre o pêndulo, sua oscilação pode ser modelada segundo a equação diferencial
não-linear (1)
g
sin (θ(t))
(1)
L
onde θ(t) é o ângulo que a massa faz com a vertical, assumindo que θ = 0 é a posição de equilíbrio. Note que
a equação diferencial dada por (1) é uma equação não-linear. Tome L = 0.3 m e g = 9.82 m/s2 , e imponha
como condição inicial θ(0) = θ0 e θ0 (0) = 0. Se o ângulo inicial θ0 está próximo de zero, então é possível usar a
seguinte aproximação sin (θ(t)) ≈ θ(t), simplificando a equação (1), obtendo a equação linearizada
θ00 (t) = −
g
θ00 (t) = − θ(t),
L
(2)
que é uma equção mais simples de se resolver.
a) Calcule o período de oscilação, T0 , do pêndulo para a equação simplificada (2);
b) Caso o valor de θ0 não seja tão pequeno, a simplificação não faz sentido, e portanto o período de oscilação
do pêndulo se torna
s
L
K(sin2 (θ0 /2)),
g
T (θ0 ) = 4
onde K(s2 ) é uma integral elíptica do primeiro tipo, dada por
Z 1
dt
2
√
√
K(s ) =
.
2
1 − s t2 1 − t2
0
Calcule numericamente uma aproximação e grafique T (θ0 ) para diferentes valores de θ0 , entre [0, 0.999π],
usando métodos de quadratura e algum outro método de integração numérica, com uma tolerância de 10−6
para o erro.
c) Verifique que, para pequenos valores de θ0 as equações (1) e (2) tem aproximadamente o mesmo período.
Estime o maior valor de θ0 tal que os períodos de oscilação tenham uma diferença relativa de no máximo
10−1 . Estude o que acontece com a razão T0 /T (θ0 ) quando θ0 → π.
d) Aproxime numericamente as soluções para o modelo não linear para diversos valores de θ0 , tanto próximos
de 0 quanto de π. Grafique seus resultados.
2
Atividade 02
Neste problema, vamos apresentar um outro método de resolução de um Problema de Valor de Contorno (PVC),
utilizando splines cúbicas por partes. Para isto, antes de introduzir o problema, apresentaremos tipo de função
que irá servir como base para o espaço das funções de classe C 2 [a, b] (funções contínuas que tem a primeira e a
segunda derivadas contínuas no intervalo dado), onde vamos querer encontrar a solução do PVC. [10]
Considere uma partição do intervalo dado por [a, b], em n subintervalos igualmente espaçados,
π : a = t0 < t1 < · · · < tn = b, tal que tı − tı−1 = h, para todo ı = 1 : n. A B-spline B (t) associada
ao ponto t ( = 1 : n) é construída impondo-se as seguintes condições
• é uma função cúbica por partes;
• tem suporte compacto, isto é, assumem valores não nulos em apenas um subintervalo. No caso das splines
cúbicas, este suporte tem comprimento 4h;
• B (t) ∈ C 2 [a, b]
Além disso, as B-splines são construídas a partir de um sistema, criado das condições impostas. Escreva de
maneira simplificada (não precisa resolver!) o sistema associado a uma B-spline [8].
Entretanto, os valores da função e de suas derivadas em alguns nós são dados na Tabela 1. Desenvolva a
expressão para uma B-spline usando os dados da tabela, e faça um gráfico da função por partes obtida.
Tabela 1: Valores para a constução das B-splines
t−2
t−1
t0
t1
t2
B(t)
0
1
4
1
0
B 0 (t)
00
B (t)
0
3/h
0
2
−3/h
0
2
−12/h
6/h
2
6/h
0
0
As funções obtidas formam o conjunto B = {B−1 (t), B0 (t), . . . , Bn+1 (t)}. Mostre que essas funções são
linearmente independentes em [a, b]. Se temos apenas n + 1 pontos, por que tivemos que adicionar mais duas
funções? Explique e faça uma representação geométrica.
Agora, o PVC considerado está mostrado em (3).


y 00 (t) + p(t)y 0 (t) + q(t)y(t) = f (t)



y(a) = α
(3)



y(b) = β,
onde as funções p(t), q(t) e f (t) são funções dadas, p(t) e q(t) são contínuas e a < b, α e β são números reais.
Para a aproximação, fazemos
n+1
X
y(t) ≈
ck Bk (t).
k=−1
Impondo as condições no PVC (3), e lembrando que a B-spline tem suporte compacto (vale zero fora do intervalo
que contém os cinco pontos que a define) encontre o sistema de equações que deve ser resolvido para encontrar
os coeficientes ck .
Por último, resolva o PVC abaixo usando o método encontrado, aproveitando ao máximo a estrutura do sistema
a ser resolvido. Altere o número de nós e grafique o resultado, comparando com a solução exata do PVC.
3

2

y 00 (t) + 4ty 0 (t) + 4π 2 y(t) = −2e−t



y(0) = 0



y(1) = e−1 ,
2
cuja solução exata é dada por y(t) = te−t cos(2πt).
4
2t3 + t cos(2πt) + 2π sin(2πt)
(4)
Atividade 03
No século XIX, o matemático A. Cayley publicou artigos sobre o método de Newton para encontrar os zeros
de uma função complexa. Ele tentou dividir o plano em regiões tal que, dada uma aproximação incial z0 , o
método convergisse para uma determinada raiz da equação considerada. Ele conseguiu determinar as regiões
para uma classe específica de funções (da forma f (z) = z 2 − c), porém não obteve sucesso para outras funções,
em decorrência dos cálculos de grande complexidade pra época. Já no início do século XX, outros dois autores se
dedicaram ao problema, G. Julia e P. Fatou, provando que, dependendo da função considerada, existem sequências
que convergem, outras que divergem ou então tem um comportamento periódico. Atualmente, é conhecido que
o método de Newton, quando aplicado a algumas funções complexas pode ter um comportamento fractal. Boas
referências para este tópico podem ser encontradas em [4, 6, 11]. Note ainda que, aplicar o método de Newton em
uma função com uma variável complexa é equivalente a aplicar o método de Newton para resolver um sistema
não linear, de duas variáveis (verifique!).
Para o conjunto de funções abaixo faça o que se pede nos itens a seguir
• Encontre a solução do sistema não linear usando o software de sua preferência;
• Escolha um retângulo no plano que contenha algumas raízes do sistema;
• Discretize esse retângulo em uma malha m × n (podendo ter m = n);
• Escreva um algoritmo em alguma linguagem de programação para resolver um sistema não linear usando
o método de Newton;
• No seu programa, associe uma cor a cada raiz, ou seja, se o chute inicial convergiu para uma raiz, pinte o
ponto (chute inicial) com a cor da raiz obtida. Associe também uma cor à falha do método;
• Rode seu programa para todos os pontos da malha, gerando uma figura;
• Se necessário, reduza o retângulo escolhido, para uma melhor visualização do fractal.
Não se esqueça dos critérios de parada necessários na implementação: Um número máximo de iterações e
também uma tolerância para se declarar que a raiz foi atingida.
(
2 − 2x + x3 − 3xy 2 = 0
1.
−2y + 3x2 y − y 3 = 0
( 5
x − 10x3 y 2 + 9x2 y + 5x2 + 5xy 4 − 4xy + 3x − 3y 3 − 5y 2 + 1 = 0
2.
5x4 y − 3x3 − 10x2 y 3 + 2x2 + 9xy 2 + 10xy + y 5 − 2y 2 + 3y = 0
(
sin (x) cosh (y) = 0
3.
cos (x) sinh (y) = 0
Crie algum sistema de equações e tente encontrar algum fractal interessante.
5
Atividade 04
A seca que atinge o estado de São Paulo é um assunto bastante explorado pela mídia ultimamente, principalmente
quando se trata dos níveis do sistema Cantareira. A intenção deste exercício é estudar um pouco o comportamento
do nível do reservatório da Cantareira nos últimos anos. Diariamente é publicado no site da SABESP [1] o nível
das represas, assim como o índice pluviométrico diário e os acumulados por mês.
O arquivo [2] disponibiliza a data (no formado AAAAMMDD). Carregue o arquivo e faça um gráfico dos dias
pelo nível do reservatório (que é dado em porcentagem). Fixe um intervalo e aplique uma transformação nos dias
para que eles se encaixem no intervalo. Qual transformação você usou?
a) Primeiramente, ajuste (via quadrados mínimos) uma reta aos dados tabelados e depois um polinômio de
grau 4. Para cada um dos casos, responda:
1. Qual o valor do resíduo obtido?
2. Este ajuste está refletindo bem o comportamento do nível da represa?
b) Proponha agora algum outra curva para o ajuste dos dados. Responda as questões acima para a curva que
você escolheu.
c) Note que, no ano de 2013 o nível da represa apresentou queda na maior parte do ano, e em 2014 foi
feito uso de uma cota do volume morto no sistema Cantareira. Selecione do arquivo os dados a partir de
01 de janeiro de 2013 e restrinja seus dados aos dias antes do uso do volume morto. Ajuste um polinômio de grau quatro a estes dados. Se não fosse possível fazer uso do volume morto, quando o sistema
ficaria com 0% de capacidade? Resolva numericamente com um método para zeros de função. Compare a
data obtida com os dados, ou seja, subtraia a porcentagem adicionada (correspondente à primeira cota do
volume morto) e encontre o dia que a capacidade zeraria.
d) Pegue os dados a partir da data em que foi acionada a segunda cota do volume morto (segundo “salto” dos
dados). Ajuste os dados por uma reta (sempre no sentido de Quadrados Mínimos) e encontre a data em que
a represa atingiria a capacidade de 0%.
6
Atividade 05
Considere uma barra homogênea de comprimento L, de forma que a largura da barra seja desprezível com relação
ao seu comprimento. Seja u(x, t) uma função que descreve a temperatura do ponto x ∈ [0, L] no instante t > 0.
Se colocarmos uma fonte de calor no ponto x = 0, como deve se comportar a temperatura na barra com o passar
do tempo? Ou então, considerando que existe uma função que descreve a temperatura da barra em um tempo
inicial, com ou sem fonte de calor, como se comporta a temperatura em relação ao tempo?
O problema proposto acima pode ser modelado como uma equação diferencial. Entretanto, temos mais de uma
variável, e assim temos um problema de uma Equação Diferencial Parcial, pois queremos avaliar o comportamento
da função u(x, t) tanto em relação ao tempo quanto ao espaço. Pode ser demonstrado que a função u(x, t) deve
ser solução da seguinte equação diferencial
∂ 2u
∂u
(x, t) = γ 2 (x, t) + ψ(x, t),
∂t
∂x
(5)
onde γ é o coeficiente de condução de calor e ψ(x, t) é uma função que representa a fonte de calor. Além disso,
devemos impor condições iniciais e de contorno sobre o problema.

u(x, 0) = f (x)



u(0, t) = p(t)



u(L, t) = q(t).
Este tipo de condição inicial dado é chamado de Condição de Dirichlet. Se, por exemplo, um dos lados da barra
estiver isolado e não conduzir calor, trocarmos a condição de contorno correspondente por uma condição do tipo
ux (x, t) = 0 (ux é a derivada de u(x, t) com relação a x), chamada de Condição de Neumann.
Uma maneira simples de se resolver problemas deste tipo é usando um esquema de diferenças finitas, fazendo
discretizações no tempo e no espaço [7]. Para isto, use
u(xi , tn ) ≈ uni .
1. Aproxime, usando diferenças finitas, as derivadas temporal e espacial do problema, utilizando:
• diferenças avançadas no tempo.
• diferenças atrasadas no tempo.
Solução:
• Para aproximar a derivada temporal usando diferenças avançadas, fazemos
uj+1
− uji
i
,
ut (xi , tj ) ≈
k
(6)
onde k é o passo associado à variável temporal, e j = 0, 1, . . . .
Agora, a derivada espacial será sempre aproximada por
uji+1 − 2uji + uji−1
uxx (xi , tj ) ≈
,
h2
onde h é tal que xı = ıh, para ı = 0, 1, . . . , m, sendo m um inteiro positivo.
Assim, podemos substituir (7) e (6) em (5) e obter
uj − 2uji + uji−1
uj+1
− uji
i
= γ i+1
+ ψij .
k
h2
7
(7)
Podemos, então, explicitar uj+1
em termo das outras variáveis, encontrando
i
uj+1
i
=
k
k
1 − 2 2 γ uji + γ 2 (uji+1 + uji−1 ) + kψij , i = 1, . . . , m, j > 0
h
h
(8)
• Para aproximar a derivada temporal usando diferenças atrasadas, fazemos
ut (xi , tj ) ≈
uji − uj−1
i
,
k
(9)
Substituindo (7) e (9) em (5), temos
uj − 2uji + uji−1
uji − uj−1
i
= γ i+1
+ ψij .
k
h2
Podemos, então, isolar uij−1 e escrevê-la em termo das outras variáveis, encontrando
uij−1
+
kψij
k j
k
k
= −γ 2 ui−1 + 1 + 2 2 γ uji − γ 2 uji+1 , i = 1, . . . , m, j > 0
h
h
h
(10)
Note que, pelo primeiro método (chamado de método explícito), conhecidas as aproximações uji para j fixo,
podemos obter facilmente as aproximações para o passo de tempo seguinte (uj+1
). Já pelo segundo método
i
(chamado de método implícito), se conhecemos uj−1
,
precisamos
resolver
um
sistema
linear para obter uji .
i
2. Escreva as condições inicial e de contorno de Dirichlet em termos da sua discretização.
Solução: A condição incial nos dá que u(x, 0) = f (x), ou seja, u0i = f (xi ), i = 0, . . . , m. Assim, dado
um valor para h = L/m, temos um vetor inicial u0 , cuja componente i é o valor uji , de dimensão m − 1.
Note ainda que, a medida que h → 0, temos que m → ∞ e portanto todos os pontos do intervalo [0, L]
tendem a ser atingidos.
Para as condições de contorno (do tipo Dirichlet), temos que uj0 = p(tj ) e ujm = q(tj ).
3. Escreva as equações obtidas na forma matricial.
Solução: Para o primeiro caso, quando temos um esquema explícito, não obtemos um sistema linear. Conseguimos representar o vetor uj+1 como a seguinte equação matricial
uj+1 = Auj + kψ j ,
onde A é a matriz

1 − 2λ
λ
0
···
 λ
1 − 2λ
λ
0

 0
λ
1
−
2λ
λ
A=
 ..
..
..
..
 .
.
.
.
0
0
···
0
(11)
···
···
···
..
.
0
0
0
..
.
λ
1 − 2λ




,


e λ = γk/h2 . Para que este método seja estável a seguinte relação deve ser satisfeita (veja em [3]):
γ
1
k
< .
2
h
2
Agora, para o segundo caso, temos que resolver um sistema linear, cuja forma matricial é dada por
8
Auj = uj−1 + kf j ,
onde A é a matriz

1 + 2λ
−λ
0
···
 −λ
1 + 2λ
−λ
0

 0
−λ
1
+
2λ
−λ
A=
 ..
..
..
..
 .
.
.
.
0
0
···
0
(12)
···
···
···
..
.
0
0
0
..
.
−λ 1 + 2λ




,


j
+λq j )> . Este segundo método tem a vantagem de ser incondicionalmente
e f j = (ψ1j +λpj , ψ2j , . . . , ψm−1
estável (veja [3]).
Aplique cada um dos métodos descritos acima para resolver numericamente o seguinte problema do calor:

∂u
∂ 2u


(x, t) =
(x, t) + 2


∂t
∂x2





u(x, 0) = sin (πx) + x(1 − x)
u(0, t) = 0





u(1, t) = 0





x ∈ [0, 1].
a) Para o método explícito, dado h = 10−2 , use dois valores para k: um tal que a condição de estabilidade seja
satisfeito e outro em que não seja.
b) Para o método implícito, use h = 10−3 e k = 10−3 , e resolva o sistema linear usando um método direto
(fatoração LU ou uma variante que aproveite a estrutura da matriz) e algum método iterativo (Gauss-Seidel,
por exemplo). Compare os resultados obtidos.
Para cada um dos itens, grafique a solução para alguns valores de tempo.
9
1.4
1.4
1.2
1.2
1.0
1.0
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.0
0.2
0.4
0.6
0.8
1.0
0.0
(a) Solução para t = 0 s
1.4
1.2
1.2
1.0
1.0
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.2
0.4
0.6
0.8
1.0
0.0
(c) Solução para t = 0.1 s
1.4
1.2
1.2
1.0
1.0
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.2
0.4
0.6
0.8
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0
(d) Solução para t = 0.2 s
1.4
0.0
0.4
(b) Solução para t = 0.03 s
1.4
0.0
0.2
1.0
0.0
(e) Solução para t = 0.5 s
0.2
0.4
0.6
0.8
(f) Solução para t = 1 s
10
1.0
Referências
[1] http://www2.sabesp.com.br/mananciais/DivulgacaoSiteSabesp.aspx. último acesso em 28-012015.
[2] http://www.ime.unicamp.br/~ms211-cursao/exercicios. último acesso em 25-02-2015.
[3] M. C. C. Cunha, Métodos Numéricos, Editora da Unicamp, 2000.
[4] A. L. S. da Cunha Freitas, Um estudo do método de newton fractal para resolução de equações, Revista da
Graduação, 4 (2011).
[5] G. Dahlquist and A. Björk, Numerical Methods, Prentice-Hall, Englewood Cliffs, NJ, 1974.
[6] L. T. dos Santos, Sistemas não lineares e fractais, Matemática Universitária, 15 (1993), pp. 102–116.
[7] R. J. LeVeque, Finite difference methods for ordinary and partial differential equations: steady-state and
time-dependent problems, Siam, 2007.
[8] P. M. Prenter, Splines and Variational Methods, A Wiley-Interscience Publication, New York, 1989.
[9] M. A. G. Ruggiero and V. L. d. R. Lopes, Cálculo numérico: aspectos teóricos e computacionais, Makron
Books do Brasil, 1997.
[10] M. H. Schultz, Spline analysis, Prentice-Hall Englewood Cliffs, NJ, 1973.
[11] S. Tatham. http://www.chiark.greenend.org.uk/~sgtatham/newton/. último acesso em 04-02-2015.
11