Notas de introdução ao Mathematica (LEFT+LEGI) 2002/03
por S. S. Valtchev ( www.math.ist.utl.pt/~ssv )

AULA 2

Programação imperativa em Mathematica (ciclos Do, While e For):

(* 
          Ciclo do tipo Do .  
              Calcula o factorial de um dado número inteiro *) 
ffDo = 
          Function[{n}, 
Module[{saida, 
                i}, 
saida = 1; 
Do[
                saida = saida*i, {i, 1, 
                  n}]; 
saida
]];

ffDo[5]

?Do

ffWhile = 
          Function[{n}, 
Module[{saida, 
                i}, 
saida = 1; 
i = 
                1; 
While[
                i ≤ n, 
saida = 
                  saida*i; 
i++;
]; 
saida
]];

ffWhile[5]

?While

ffFor = 
          Function[{n}, 
Module[{saida, 
                i}, 
saida = 1; 
For[
                i = 1, i ≤ 
                  n, i++, saida = 
                    saida*i;]; 
saida
]];

ffFor[5]

?For

Timing[ffDo[
            12345];]
Timing[ffWhile[
          12345];]
Timing[ffFor[
          12345];]

Apply[Plus, DigitCount[12345!]]  (* 
          o número 12345! tem 45151 digitos!!!! *)

Exercícios:

Ex. 1)	Implemente em Mathematica o Método do Ponto Fixo através duma função mpf[ g, x0, eps , nmax ] com 4 argumentos, onde g é a função iteradora do método, x0 é a iterada inicial para o método e eps é a tolerância, i.e. o método calcula as sucessivas iteradas (Subscript[x, n + 1]) = g((Subscript[x, n])) até que a condição |(Subscript[x, n + 1]) - 
          (Subscript[x, n])|<eps for verificada ou até um número máximo de iterações nmax. Se a condição de tolerância for verificada antes do número máximo de iterações a função deve retornar como saida a aproximação do ponto fixo com 16 casas decimais  e o número de iteradas, caso contrário, uma mensagem de erro deve ser apresentada. 
	Considere como exemplo a seguinte função iteradora g(x) = Cos[x]. Relembrando que os pontos fixos de uma função são dados pela intersecção da função g com a bissectriz do 1º e 3º quadrantes, determine um intervalo de comprimento 0.5 que contêm o único ponto fixo da função dada. Aplique a função que definiu para obter uma aproximação da raiz, considerando eps=10E-10, x0 pertencente ao intervalo que determinou e nmax=100. Utilisando o comando FindRoot do Mathematica determine o ponto fixo com uma precisao de 16 (use a opção AccuracyGoal do comando referido). Compare os valores obtidos pelo metodo implementado e pelo FindRoot em termos de erro relativo.

mpf=Function[{g, x0, eps, nmax},
Module[{a, b, n},
a = x0;   
b = g[a];    
n = 1;    
While[Abs[a - b] ≥ eps && n < nmax,    
a = b;

                      b = g[a];

                      n++;    
];
If[
                    n == 
                      nmax, 
Print["Foi atingido o n\
úmero máximo permitido de \
iterações!"];, 
Print["Depois de \
n=", n, " iteradas, foi obtida a seguinte aproximação \
do ponto fixo z=", SetPrecision[b, 16]];];
]];

g = Function[x, Cos[x]];

Plot[{x, g[x]}, {x, 0, Pi/2}, Frame -> True];

(* o intervalo pedido pode ser escolhido como [0.5,  1.0] *)

mpf[g, 0.4, 10^-10, 100]

z=x/.SetPrecision[
          FindRoot[g[x] == x, {x, 0.4}, AccuracyGoal -> 16], 16] (* valor do ponto fixo com 16 casas decimais obtido através da rotina FindRoot[] do Mathematica  *)

Abs[(z - 0.7390851331840241)/
          z]

Ex. 2)	Altere a função anterior ( mpf[g,x0,eps,nmax] ) de forma a guardar uma lista de pares {aproximação, erro absoluto}. Os erros absolutos devem ser exactos, para isso use o comando FindRoot[] do Mathematica, especificando a precisão necessária. Trace o gráfico do erro em função ao índice da iterada. Evidencie a convergência do método, tracando a órbita do ponto fixo.

mpfnew=Function[{g, x0, eps, nmax},
Module[{a, b, n},(* a lista dos erros e a solução sol não aparecem nas variáveis locais! O seu valor vai ser guardado para poder ser pedido fora da função *)
sol = 
                    x /. FindRoot[g[x] == x, {x, x0}, 
                        AccuracyGoal -> 16];
erros = {{x0, Abs[x0 - sol]}};  
a = x0;
b = g[a];
n = 1;
AppendTo[
                    erros, {b, Abs[b - sol]}];  
While[
                    Abs[a - b] ≥ eps && 
                      n < nmax, 
a = 
                      b; 
b = 
                      g[a]; 
AppendTo[
                      erros, {b, 
                        Abs[b - 
                            sol]}]; 
n++; 
];
If[
                    n == 
                      nmax, 
Print["Foi atingido o n\
úmero máximo permitido de \
iterações!"];, 
Print["Depois de \
n=", n, " iteradas, foi obtida a seguinte aproximação \
do ponto fixo z=", SetPrecision[b, 16]];];
]];

g = Function[x, Cos[x]];   
x0 = 0.1;   
mpfnew[g, x0, 10^-5, 35]

ListPlot[erros, PlotJoined -> True, PlotRange -> All, 
        Frame -> True, GridLines -> {{sol}, {}}, 
        Axes -> False]

aux = 
          Table[{i - 1, erros[[i, 2]]}, {i, 1, Length[erros]}];  
ListPlot[
      aux, PlotJoined -> True, PlotRange -> All, 
      Frame -> True, Axes -> False]

orbita = {{x0, 
              0}};
Do[
AppendTo[
          orbita, {erros[[i - 1, 1]], 
            erros[[i, 1]]}]; 
AppendTo[
          orbita, {erros[[i, 1]], 
            erros[[i, 1]]}];
, {i, 2, 
          Length[erros]}];
graf1 = 
        Plot[{x, g[x]}, {x, 0, Pi/2}, 
          PlotStyle -> {RGBColor[0, 0, 0], RGBColor[1, 0, 0]}, 
          DisplayFunction -> 
            Identity];
graf2 = 
        ListPlot[orbita, PlotJoined -> True, 
          DisplayFunction -> Identity];
Show[
        graf1, graf2, DisplayFunction -> $DisplayFunction, 
        Frame -> True, AspectRatio -> Automatic];

Table[
        erros[[i, 2]]/erros[[i - 1, 2]], {i, 20, Length[erros]}] (* cálculo de uma aproximação do coeficiente assimptotico de convergência do método do ponto fixo de ordem p=1.  Definição-- >  
                (Subscript[K, ∞])=\(\(|\)\(e\_\(n + 
                      1\)\)\(|\)\(\(/\)\(\(|\)\(e\_n\)\( | \^p\)\)\)\)  *)

Ex. 3)	Implemente em Mathematica uma função mat[ func, dim , intx, inty] que constrói uma matriz quadrada de dimensão dim. Cada componente ij da matriz é dada pelo valor da função func no ponto {xi,yj}, onde xi e yj são pontos igualmente espaçados no intervalo intx=[ax,bx] (argumento da função mat) e inty=[ay,by] (argumento da função mat), respectivamente.

mat=Function[{func, dim, intx, inty},
Module[{hx, hy, aux, x, y},
hx = 
                    N[(intx[[2]] - intx[[1]])/(dim - 1)];   
hy = 
                    N[(inty[[2]] - inty[[1]])/(dim - 1)];    
aux = Table[0.0, {dim}, {dim}];    
Do[
Do[
x = 
                              intx[[1]] + (j - 1)*
                                  hx;

                          y = 
                              inty[[2]] - (i - 1)*
                                  hy; 

                          aux[[i, j]] = func[x, y];  (* 
                              preencher a matriz entrada por entrada *) 
,{j, 1, dim}]
,{i, 1, dim}];
aux   
]];

f[x_, y_] := 
          Exp[-(x^2 + y^2)]*
            Sin[(x^3 + y^3)*2*
                Pi];
intx = {-Pi/4, 
          Pi/4};
inty = {-Pi/4, 
          Pi/4};
npontos = 50; 
Plot3D[
        f[x, y], {x, intx[[1]], intx[[2]]}, {y, inty[[1]], 
          inty[[2]]}, PlotPoints -> npontos, Mesh -> False, 
        PlotRange -> All];
DensityPlot[
        f[x, y], {x, intx[[1]], intx[[2]]}, {y, inty[[1]], 
          inty[[2]]}, PlotPoints -> npontos, Mesh -> False, 
        ColorFunction -> Hue];

ppp = mat[f, npontos, intx, inty];   
ppp = 
        Reverse[ppp];
ListPlot3D[ppp, 
      PlotRange -> All, Mesh -> False]

Dimensions[ppp]

Det[ppp]

Ex. 4)	Implemente em Mathematica uma função que calcula o número de condição de uma matriz, na norma infinito (\(TraditionalForm\`\(cond\_ ∝ \) = \(\(||\)\(A\)\( || \_ ∝ \)\(||\)\(A\^\(-1\)\)\( || \_ ∝ \)\)\)). Se a matriz introduzida não for quadrada a função deve devolver uma mensagem de erro. Para começar implemente uma função chamada norma[] que calcula a norma infinito duma matriz introduzida pelo utilizador. Exemplifique, utilizando a matriz do exercício anterior e uma outra matriz à sua escolha.

norma=Function[{matriz},
Module[{dim, aux, auxaux},
dim = Dimensions[matriz];
aux = Abs[matriz];   
auxaux = 
                    Table[Apply[Plus, aux[[i]]], {i, 1, dim[[1]]}]; (* a  s somas por linhas são guardadas no vector auxaux *)
Max[auxaux]
]];

?Apply

norma[ppp]

cond = 
        Function[{matriz}, 
Module[{inv, 
              dim}, 
dim = 
              Dimensions[matriz]; 
If[
              dim[[1]] == dim[[2]], 
inv = 
                Inverse[matriz]; 
norma[matriz]*
                norma[inv], 
Print["A matriz nao \
é quadrada!"]]
]];

cond[ppp]

(* O númeor de condição da matriz ppp é muito grande o que implica que a resolução dum sistem linear com esta matriz vai ser um problema muito mal condicionado *)

(* consideramos por exemplo uma matriz A tridiagonal (dim = 10) com valor 1 no diagonal principal e com 0  .1 nos dois diagonais adjacentes *)
dim = 15;
A = 
        IdentityMatrix[dim];
Do[
        Do[If[j == i + 1 || j == i - 1, A[[i, j]] = 0.1], {j, 
            1, dim}], {i, 1, dim}];
MatrixForm[
      A]

cond[A]