Uma Introdução ao Maple V - joinville.udesc.br Uma Introdução ao Maple V 2 Álgebra de Matrizes e Vetores A maioria da álgebra linear está incluída no pacote de rotinas do Maple

Uma Introdução ao Maple Vpor Jun S. O. Fonseca.

SumárioIniciando ........................................................................................................ 1Álgebra de Matrizes e Vetores ...................................................................... 2

Definir uma matriz A 3 por 4 ............................................................ 2Encontrando os autovalores de uma matriz quadrada AAT.............. 2Matrizes Complexas .......................................................................... 3Álgebra de vetores............................................................................. 4Álgebra de matrizes ........................................................................... 5

Solução de Equações Algébricas................................................................... 6Funções e sua Representação Gráfica ........................................................... 7

Gráficos ............................................................................................. 7Encontrando Raízes e Mínimos......................................................... 9

Comandos de Cálculo e de Programação ...................................................... 9Integração: ......................................................................................... 9Diferenciação e Integração:............................................................... 10Séries de Fourier:............................................................................... 11Ajuste de Curvas................................................................................ 11

Exemplos de Mecânica dos Sólidos I ............................................................ 13Reações de Apoio e Deflexão de uma Viga com Carga Distribuída. 14Outro problema de deflexão de uma viga.......................................... 15

Exemplos de Análise Estrutural Avançada ................................................... 18Applicação do Teorema de Castigliano............................................. 18Uma Aplicação do Método de Rayleigh-Ritz ................................... 18Segunda Aplicação do Método de Rayleigh-Ritz ............................. 20Exemplo de um Pórtico Plano ........................................................... 22

Exemplos em Análise Estrutural Avançada .................................................. 23Exemplo de Flexão: ........................................................................... 25

Iniciando

Para rodar o Maple em um computador em que este programa esteja instalado, deve-se:• procurar seu ícone nos menus ou grupos de programas (MS Windows, SGI Irix, HP e

SUN CDE, Apple MacOS, Linux Gnome ou KDE) e clicá-lo;• digitar xmaple no “prompt” em um “shell” (na maioria dos sistemas Unix, Linux, e

FreeBSD)• Digitar maple no “prompt” do DOS prompt, se a via (“path”) estiver certo.

Uma vez que o Maple estiver rodando, digita-se os comandos no prompt >. Umcomando pode ocupar muitas linhas; o final do comando é reconhecido pelo caracter definalização “;” ou “:”. A finalização com dois pontos faz com que o Maple não imprima asaída do comando. Os resultados dos últimos comandos são armazenados na memória epodem ser recuperados com os símbolos %, %%, e %%%.

O comando mais útil para os iniciantes é o de ajuda. A ajuda pode ser acionada atravésdos menus ou digitando-se o caracter ?, que não precisa finalização. Uma boa dica é utilizaros tutoriais. Para conseguir ajuda sobre um comando específico, como por exemplo solvepode-se digitar ?solve ou selecionar ajuda com o botão direito do mouse.

Uma Introdução ao Maple V

2

Álgebra de Matrizes e Vetores

A maioria da álgebra linear está incluída no pacote de rotinas do Maple chamadolinalg. Desta maneira, antes de começar lidar com matrizes e vetores, deve-se carregar estepacote com o comando:

> with(linalg):

Este é uma linha de comando típica do Maple. O comando em si é o with, que é seguido porseus argumentos entre parênteses (linalg), e por um caracter de finalização. Este caracter podeser tanto dois pontos “:” ou ponto e vírgula “;”, conforme se queira que o Maple imprima ounão o resultado do comando.

Definir uma matriz A 3 por 4Pode-se definir um conjunto de números ou símbolos em uma matriz 3x4:

> A:=matrix(3,4,[1,2,3,4,5,4,7,8,-4,-3,-2,-1]); [ 1 2 3 4 ] [ ] A := [ 5 4 7 8 ] [ ] [ -4 -3 -2 -1 ]

Para imprimir o valor de uma variável, pode-se usar tanto o comando print> print(A);

[ 1 2 3 4 ] [ ] [ 5 4 7 8 ] [ ] [ -4 -3 -2 -1 ]

como um dos comandos de avaliação eval. Neste caso, usa-se a avaliação de matrizes evalm.> evalm(A);

[ 1 2 3 4 ] [ ] [ 5 4 7 8 ] [ ] [ -4 -3 -2 -1 ]

Os comandos de avaliação são muito importantes no Maple. Existem muitos comandosinertes, isto é, comandos que não são executados imediatamente, mas somente quando se pedeao Maple uma avaliação da expressão. O leitor notará o uso dos comandos inertes ao longo dotexto. Estes comandos normalmente iniciam com letra maiúscula.

Encontrando os autovalores de uma matriz quadrada AAT

Calcula-se agora os autovalores de uma matriz quadrada obtida pela multiplicaçãoAAT. As maneiras de se fazer isto no Maple são:

> V:=eigenvals(multiply(A,transpose(A)));

ou> V:=eigenvals(A &* transpose(A));

1/2 V := %3 + %2 + 214/3, - 1/2 %3 - 1/2 %2 + 214/3 + 1/2 3 (%3 - %2)

1/2 - 1/2 %3 - 1/2 %2 + 214/3 - 1/2 3 (%3 - %2) I

1/2 1/2%1 := I 19226186 3

/7442560 \1/3%2 := |------- - 16/3 %1| \ 27 /

/7442560 \1/3%3 := |------- + 16/3 %1| \ 27 /


3

Esta é uma computação simbólica, imune a erros de arredondamento. No entanto, estetipo de computação requer muitas operações (e muito tempo) de computador. Na maioria dasvezes, é mais prático calcular os autovalores numéricamente, especialmente se a matriz formaior que 4x4. Uma possibilidade é converter a matriz para uma matriz de números de“ponto-flutuante” ao invés de números exatos usando o comando de avaliação evalf, que faráeigenvals calcular valores de ponto-flutuante:

> k:=eigenvals(evalf(evalm(A &* transpose(A))));

k := 201.8047330, .598203057, 11.59706392

mas isto toma muito tempo. A melhor maneira é calcular numéricamente desde o início,usando o método QR implementado no comando Eigenvals do Maple. Note que estecomando é inerte; o cálculo acontecerá somente com o comando de avaliação.

> k:=evalf(Eigenvals(multiply(A,transpose(A))));

k := [ 201.8047330, .598203057, 11.59706392 ]

Os autovetores de uma matriz também podem ser obtidos pelo Maple, através do comando> v:=eigenvects(multiply(A,transpose(A)));

v :=

2 469 4584 2 397 7325 [%1, 1, {[ - 9/5356 %1 + ---- %1 - ----, 5/8034 %1 - ---- %1 + ----, 1 ]}] 1339 1339 2678 4017

2 3%1 := RootOf(- 1400 + 2468 _Z - 214 _Z + _Z )

ou numericamente através do comando inerte Eigenvals, que tem um argumento adicionalque é onde os autovetores serão armazenados:

> k:=evalf(Eigenvals(multiply(A,transpose(A))),vecs));

k := [ 201.8047334, .59820303, 11.59706393 ]

> evalm(vecs);

+ .3917147767 .8136191924 -.3722840080 + ¦ ¦ ¦ .9183834076 -.4391516693 -.1696430687 ¦ ¦ ¦ + -.3342573782 -.2531071924 -.9023780667 +

Uma matriz chamada vecs (ou qualquer outro nome escolhido) será criada tendo osautovetores como colunas. Estas operações são úteis, por exemplo, para o cálculo de tensões edireções principais.

Matrizes ComplexasO Maple define o símbolo I como um sinônimo da raiz de -1 (sqrt(-1)). Pode-se utilizar estesímbolo para definir uma matriz complexa, através, por exemplo, de duas matrizes:

> A:=matrix(3,3,[1,2,3,3,2,1,-1,1,-1]); [ 1 2 3 ] [ 3 2 1 ] [ -1 1 -1 ]

> B:=matrix(3,3,[1,0,1,-1,-1,1,4,1,2]); [ 1 0 1 ] [ -1 -1 1 ] [ 4 1 2 ]

> AB:=matadd(A,I*B);

or> AB:=evalm(A+I*B);

[ 1 + I 2 3 + I ] [ 3 - I 2 - I 1 + I ] [ -1 + 4 I 1 + I -1 + 2 I ]

Pode-se também utilizar diretamente os números complexos:


4

> AB:=matrix(3,3,[1+I,2,3+I,3-I,2-I,1+I,-1+4*I,1+I,-1+2*I]); [ 1 + I 2 3 + I ] [ 3 - I 2 - I 1 + I ] [ -1 + 4 I 1 + I -1 + 2 I ]

Para se extrair a parte real e imaginária de uma matriz complexa, pode-se utilizar o comandomap, que aplica uma operação em todos os termos de uma matriz, e o comando de avaliaçãode variáveis complexas evalc, e a definição de uma função ”->”.

> map(x -> evalc(Re(x)),AB); [ 1 2 3 ] [ 3 2 1 ] [ -1 1 -1 ]

> map(x -> evalc(Im(x)),AB); [ 1 0 1 ] [ -1 -1 1 ] [ 4 1 2 ]

Álgebra de vetoresExiste um conjunto de operações no pacote de rotinas “linalg” que podem ser úteis para amanipulação de vetores. Tem-se, por exemplo, o comando crossprod para produto externo noR3 e dotprod e innerprod para produtos internos, sendo que este último aceita matrizescomo argumento. Cria-se inicialmente dois vetores:

> v:=vector(3); v := array(1 .. 3, [])

> w:=vector(3); w := array(1 ..3, [])

> evalm(v); [v[1], v[2], v[3])

É interessante ressaltar que o comando matrix sempre cria um arranjo bidimensional e ocomando vector sempre um arranjo unidimensional.

> evalm(dotprod(v,w)); v[1] w[1] + v[2] w[2] + v[3] w[3]

> evalm(crossprod(v,w); [ v[2] w[3] - v[3] w[2], v[3] w[1] - v[1] w[3], v[1] w[2] - v[2] w[1] ]

> evalm(innerprod(v,matrix(3,3,[1,0,2,0,3,0,2,0,2]),w)); v[1] w[1] + 2 v[3] w[1] + 3 v[2] w[2] + 2 v[1] w[3] + 2 v[3] w[3]

Algumas operações termo a termo não estão implementadas diretamente no Maple, edevem ser programadas pelo usuário através de comandos de laços ou mapeamento. Porexemplo, se for necessário dividir cada termo de um vetor pelo termo correspondente deoutro, pode-se escrever

> u:=vector(3);> for i from 1 to 3 do> u[i]:=v[i]/w[i]> od:> evalm(u);

u := array(1 .. 3, [])

v[1] v[2] v[3] [ ----, ----, ---- ] w[1] w[2] w[3]

Pode-se, alternativamente, implementar-se esta operação diretamente na definição dovetor, operando-se sobre o índice:

> t:=vector(3,i->v[i]/w[i]); v[1] v[2] v[3] t := [ ----, ----, ---- ] w[1] w[2] w[3]


5

Álgebra de matrizesO Maple usa símbolos diferentes para operações com números e com matrizes.

Entretanto, alguns dos operadores normais funcionam dentro do comando evalm, com aexcessão do operador de multiplicação * , reservado para a multiplicação por escalar. Amultiplicação de matrizes não-comutativa é representada pelo operador “&*”.

> A:=matrix(3,3):> evalm(A);

[ A[1, 1] A[1, 2] A[1, 3] ] [ ] [ A[2, 1] A[2, 2] A[2, 3] ] [ ] [ A[3, 1] A[3, 2] A[3, 3] ]

> B:=matrix(3,3):> evalm(B);

[ B[1, 1] B[1, 2] B[1, 3] ] [ ] [ B[2, 1] B[2, 2] B[2, 3] ] [ ] [ B[3, 1] B[3, 2] B[3, 3] ]

> matadd(A,B); or > evalm(A + B); [ A[1, 1] + B[1, 1] A[1, 2] + B[1, 2] A[1, 3] + B[1, 3] ] [ ] [ A[2, 1] + B[2, 1] A[2, 2] + B[2, 2] A[2, 3] + B[2, 3] ] [ ] [ A[3, 1] + B[3, 1] A[3, 2] + B[3, 2] A[3, 3] + B[3, 3] ]

> multiply(A,B); or > evalm(A &* B); [A[1, 1] B[1, 1] + A[1, 2] B[2, 1] + A[1, 3] B[3, 1],

A[1, 1] B[1, 2] + A[1, 2] B[2, 2] + A[1, 3] B[3, 2],

A[1, 1] B[1, 3] + A[1, 2] B[2, 3] + A[1, 3] B[3, 3]]

[A[2, 1] B[1, 1] + A[2, 2] B[2, 1] + A[2, 3] B[3, 1],

A[2, 1] B[1, 2] + A[2, 2] B[2, 2] + A[2, 3] B[3, 2],

A[2, 1] B[1, 3] + A[2, 2] B[2, 3] + A[2, 3] B[3, 3]]

[A[3, 1] B[1, 1] + A[3, 2] B[2, 1] + A[3, 3] B[3, 1],

A[3, 1] B[1, 2] + A[3, 2] B[2, 2] + A[3, 3] B[3, 2],

A[3, 1] B[1, 3] + A[3, 2] B[2, 3] + A[3, 3] B[3, 3]]

> multiply(A,inverse(B)); or > evalm(A / B);[(A[1, 1] B[2, 2] B[3, 3] - A[1, 1] B[2, 3] B[3, 2] - A[1, 2] B[2, 1] B[3, 3]

+ A[1, 2] B[2, 3] B[3, 1] + A[1, 3] B[2, 1] B[3, 2] - A[1, 3] B[2, 2] B[3, 1]

)/%1, - (A[1, 1] B[1, 2] B[3, 3] - A[1, 1] B[1, 3] B[3, 2]

- A[1, 2] B[1, 1] B[3, 3] + A[1, 2] B[1, 3] B[3, 1] + A[1, 3] B[1, 1] B[3, 2]

- A[1, 3] B[1, 2] B[3, 1])/%1, (A[1, 1] B[1, 2] B[2, 3]

- A[1, 1] B[1, 3] B[2, 2] - A[1, 2] B[1, 1] B[2, 3] + A[1, 2] B[1, 3] B[2, 1]

+ A[1, 3] B[1, 1] B[2, 2] - A[1, 3] B[1, 2] B[2, 1])/%1]

[(A[2, 1] B[2, 2] B[3, 3] - A[2, 1] B[2, 3] B[3, 2] - A[2, 2] B[2, 1] B[3, 3]

+ A[2, 2] B[2, 3] B[3, 1] + A[2, 3] B[2, 1] B[3, 2] - A[2, 3] B[2, 2] B[3, 1]

)/%1, - (A[2, 1] B[1, 2] B[3, 3] - A[2, 1] B[1, 3] B[3, 2]

- A[2, 2] B[1, 1] B[3, 3] + A[2, 2] B[1, 3] B[3, 1] + A[2, 3] B[1, 1] B[3, 2]

- A[2, 3] B[1, 2] B[3, 1])/%1, (A[2, 1] B[1, 2] B[2, 3]

- A[2, 1] B[1, 3] B[2, 2] - A[2, 2] B[1, 1] B[2, 3] + A[2, 2] B[1, 3] B[2, 1]

+ A[2, 3] B[1, 1] B[2, 2] - A[2, 3] B[1, 2] B[2, 1])/%1]


6

[(A[3, 1] B[2, 2] B[3, 3] - A[3, 1] B[2, 3] B[3, 2] - A[3, 2] B[2, 1] B[3, 3]

+ A[3, 2] B[2, 3] B[3, 1] + A[3, 3] B[2, 1] B[3, 2] - A[3, 3] B[2, 2] B[3, 1]

)/%1, - (A[3, 1] B[1, 2] B[3, 3] - A[3, 1] B[1, 3] B[3, 2]

- A[3, 2] B[1, 1] B[3, 3] + A[3, 2] B[1, 3] B[3, 1] + A[3, 3] B[1, 1] B[3, 2]

- A[3, 3] B[1, 2] B[3, 1])/%1, (A[3, 1] B[1, 2] B[2, 3]

- A[3, 1] B[1, 3] B[2, 2] - A[3, 2] B[1, 1] B[2, 3] + A[3, 2] B[1, 3] B[2, 1]

+ A[3, 3] B[1, 1] B[2, 2] - A[3, 3] B[1, 2] B[2, 1])/%1]

%1 := B[1, 1] B[2, 2] B[3, 3] - B[1, 1] B[2, 3] B[3, 2] - B[2, 1] B[1, 2] B[3, 3]

+ B[2, 1] B[1, 3] B[3, 2] + B[3, 1] B[1, 2] B[2, 3] - B[3, 1] B[1, 3] B[2, 2]

Similarmente à algebra de vetores, muitas operações termo a termo em matrizes tem que serdefinidas pelo usuário. Por exemplo, uma divisão termo a termo de uma matriz por outrapoderia que ser definida por

> C:=matrix(3,3,(i,j)->A[i,j]/B[i,j]); [ A[1, 1] A[1, 2] A[1, 3] ] [ ------- ------- ------- ] [ B[1, 1] B[1, 2] B[1, 3] ] [ ] [ A[2, 1] A[2, 2] A[2, 3] ] C := [ ------- ------- ------- ] [ B[2, 1] B[2, 2] B[2, 3] ] [ ] [ A[3, 1] A[3, 2] A[3, 3] ] [ ------- ------- ------- ] [ B[3, 1] B[3, 2] B[3, 3] ]

Pode-se extrair o numerador e o denominador dos termos desta matriz usando o comando mappara mapear os comandos numer e denom para cada termo:

> map(x->numer(x),C); [ A[1, 1] A[1, 2] A[1, 3] ] [ ] [ A[2, 1] A[2, 2] A[2, 3] ] [ ] [ A[3, 1] A[3, 2] A[3, 3] ]

> map(x->denom(x),C); [ B[1, 1] B[1, 2] B[1, 3] ] [ ] [ B[2, 1] B[2, 2] B[2, 3] ] [ ] [ B[3, 1] B[3, 2] B[3, 3] ]

Solução de Equações Algébricas

Pode-se resolver um sistema de equações simultâneas lineares (e não-lineares em algunscasos) simbolicamente no Maple utilizando o comando solve. Define-se as equações (ouconjunto de equações) e as variáveis a serem resolvidas:

> solve({a11*x1+a12*x2=b1, a21*x1+a22*x2=b2},{x1,x2});

a12 b2 - b1 a22 a11 b2 - a21 b1 {x1 = - -----------------, x2 = -----------------} a11 a22 - a21 a12 a11 a22 - a21 a12

Se o Maple conseguir encontrar a solução, ele devolve um conjunto de identidades que podemser aplicadas às variáveis através do comando assign, ou substituída em alguma expressãocom o comando subs, por exemplo

> a:=%:> f:= subs(a,2*x1+x2);

a12 b2 - b1 a22 a11 b2 - a21 b1 f := - 2 ----------------- + ----------------- a11 a22 - a21 a12 a11 a22 - a21 a12


7

Funções e sua Representação Gráfica

Funções podem ser definidas de duas maneiras no Maple. Para funções simples, pode serconveniente utilizar o operador “seta” ->, para uma definição direta,

> f:=x -> x-2*x+x^2+x*cos(exp(-x)); 2 f := x -> - x + x + x cos(exp(- x))

> g:=y -> y+y^2*sin(1/(1+y^2)); 2 1 g := y -> y + y sin(------) 2 1 + y

> h:=(x,y) -> f(x*y)+g(x+y); h := (x,y) -> f(x y) + g(x + y)

Funções mais complicadas podem ser definidas utilizando o comando proc como a seguir,> f1:=proc(x) if x>3 then x^2 else if x <= 3 then x-5 fi fi end;

f1 := proc(x) if 3 < x then x^2 else if x <= 3 then x-5 fi fi end

Há muitos comandos para manipulação de expressões. Por exemplo, pode-se reescrever umafunção de muitas variávies organizando-as pelas potências de x ou y,

> collect(h(x,y),x); / 2 1 \ 2 |y + sin(------------)| x | 2 | \ 1 + (x + y) /

/ 1 \ + |y cos(exp(- x y)) + 1 - y + 2 y sin(------------)| x + y | 2 | \ 1 + (x + y) /

2 1 + y sin(------------) 2 1 + (x + y)

GráficosMaple tem vários comandos implementados para traçar gráficos, e muitas opções paracontrolar a aparência. Os mais simples são plot e plot3d. Para um gráfico simples em duasdimensões, pode-se digitar

> plot(f,-4..4);> plot(g,-4..4);

Pode-se traçar gráficos paramétricos, tais como


8

> plot([f(t),g(t),t=-4..4]);

Os gráficos tridimensionais são feitos de forma similar. É possível especificar muitas dasopções dos comandos gráficos diretamente no comando plot3d, mas normalmente é mais fácilmudar os parâmetros de visualização diretamente nos menus da janela gráfica. Por exemplo,pode-se mudar as opções posição do observador e ligar as curvas de nível para fazer osegundo gráfico abaixo.

> fg:=(p,q)->f(p*q)*g(1+p-q);> plot3d(fg,-2..2,-2..2);

Existem centenas de outras opções de plotagem no pacote de rotinas plots. Estescomandos adicionais permitem a geração de uma rica gama de figuras, tais quais asespetaculares figuras de capas de livros de matemática.


9

Encontrando Raízes e MínimosÉ possível encontrar raízes de uma equação não-linear definida por uma função utilizando ocomando solve (que busca uma solução simbólica e normalmente não encontra) ou fsolve(que busca uma solução numérica e normalmente só encontra se houver um bom palpite dalocalização). Por exemplo, se o objetivo for encontrar na função abaixo um valor de x para oqual a função valha 5, faz-se

> f:=x->x-2*x+x^2+x*cos(exp(-x)):> fsolve(f(x)=5,x,1..3);

2.238906646

> fsolve(f(x)=5,x,-3..1); -1.994108139

Para se encontrar extremos (máximos e mínimos) de uma função, utiliza-se o método quedeveria ter sido aprendido em Cálculo. Normalment uma combinação de plots e fsolvesfunciona melhor para localizar o extremo.> b:=fsolve(diff(f(x),x),x,-4.1..-4); evalf(f(b)); b := -4.034475437

16.27951972

Comandos de Cálculo e de Programação

Depois de definir uma função qualquer, pode-se imprimí-la em uma variedade deformatos. O Maple tem muitas opções, desde o horrível modo texto usado neste documentoaté as formas visíveis na tela do computador ou no formato PostScript. Abaixo se usou asformas inerte Int e Diff para que o Maple imprima sem avaliar.

> f:=x->x-2*x+x^2+cos(exp(-x)); Int(f(x),x=0..x2); Diff(y,x);

É possível também converter as expressões do Maple para as linguagens FORTRAN e C, parauma possível inclusão em um programa de computador, ou na linguagem LaTEX paraprocessadores de textos científicos. Para direcionar a saída de um comando para um arquivoao invés da tela, algo bastante útil nos casos acima, utiliza-se os comandos write e writeto.

> f:=x->x-2*x+x^2+cos(exp(-x)); 2 f := x -> - x + x + cos(exp(- x))

> fortran (f(x)); readlib(C): C(f(x)); t0 = -x+x**2+cos(exp(-x)) t0 = -x+x*x+cos(exp(-x));

> latex(f(x));-x+x^{2}+\cos({e^{-x}})

Integração:Pode-se usar a integração analítica ou numérica. O operador de integração tem duas

formas diferentes, a forma ativa int, e a inerte Int. Para se utilizar uma integraçãocompletamente númerica, deve-se mandar o Maple avaliar em ponto flutuante (evalf) a formainerte da integral.


10

Existem vários comandos de controle de execução no Maple, para auxiliar nacomputação de procedimentos repetitivos. Laços de execução, por exemplo, são definidoscom os comandos for ... from ... to ... by ... do { } od: .

Apresenta-se agora uma série de comandos de séries que fizeram a fama do Maple. Ocomando taylor gera uma série de Taylor's com uma precisão especificada. Note que umasérie, com seu termo de ordem, é uma entidade completamente diferente de um polinômio.Deve-se convertê-la para um polinômio para manipulação e avaliação.

> difference:=array(1..3); difference := array(1 .. 3, [])

> exact:=evalf(Int(f(x),x=-2..2)); exact := 6.870394434

> for i from 1 to 3 do> FI:=taylor(f(x),x=0,2*i):> print(`Taylor s Series `,2*i,` = `,FI);> Poly:=convert(FI,polynom):> difference[i]:=evalf(exact - int(Poly,x=-2..2));> od:

2 Taylor s Series , 2, = , cos(1) + (- 1 + sin(1)) x + O(x )

Taylor s Series , 4, = ,

2 3 cos(1) + (- 1 + sin(1)) x + (1 - 1/2 cos(1) - 1/2 sin(1)) x + 1/2 cos(1) x

4 + O(x )

Taylor s Series , 6, = ,

2 3 cos(1) + (- 1 + sin(1)) x + (1 - 1/2 cos(1) - 1/2 sin(1)) x + 1/2 cos(1) x

4 / 23 \ 5 6 + (5/24 sin(1) - 1/4 cos(1)) x + |- --- sin(1) + 1/24 cos(1)| x + O(x ) \ 120 /

> evalm(difference); [ 4.709185210, 3.060580653, 2.545625406 ]

Diferenciação e Integração:O Maple pode realizar integração analítica se for preciso; mas deve-se considerar que

esta operação é bastante demorada. Podemos comparar a diferença entre a integração analíticae numérica. Estes resultados podem ser melhorados aumentando a precisão do comando evalf.

> diff(f(x),x);> int(f(x),x=0..1);> fiexact:=%;> fiquadrat:=evalf(Int(f(x),x=0..1));> fiquad20:=evalf(Int(f(x),x=0..1),20);> evalf(fiexact-fiquadrat);evalf(fiexact-fiquad20,20);

- 1 + 2 x + sin(exp(- x)) exp(- x)

- 1/6 - Ci(exp(-1)) + Ci(1)

fiexact := - 1/6 - Ci(exp(-1)) + Ci(1)

fiquadrat := .6271651970

fiquad20 := .62716519698330620506

-9 -.1*10 -19 -.2*10


11

Séries de Fourier:Pode-se calcular os coeficientes de uma expansão em série de Fourier através da definiçãobásica. Para isto, utiliza-se a capacidade de integração do Maple. Para fins de comparação,calculamos agora a expansão em série de Fourier da função de Heaviside (degrau unitário),que é pré-pogramada no Maple.

> f:=t->Heaviside(t); f := Heaviside

> fsn:=(int(f(x),x=-Pi..Pi))/(2*Pi)+sum(int(f(x)*sin((2*k-1)*x),x=-Pi..Pi)*sin((2*k-1)*x)/Pi,k=1..5);

sin(x) sin(3 x) sin(5 x) sin(7 x) sin(9 x) fsn := 1/2 + 2 ------ + 2/3 -------- + 2/5 -------- + 2/7 -------- + 2/9 -------- Pi Pi Pi Pi Pi

> plot(fsn,x=-Pi..Pi,title=`Fourier series`);

Poderia-se alternativamente ter definido a função degrau por proc(x) if x<0 then f:=x-> 0 else if x>=0 then f:=x->1 fi fi end:.

Ajuste de CurvasUsando o Maple, pode-se fazer ajuste de curvas sobre um conjunto de dados discretos.

Pode-se por exemplo simular um conjunto de 100 pontos usando a avaliação de funções edepois realizar um ajuste utilizando uma base de funções conhecidas. Os coeficientes de cadatermo da base podem ser encontrados por mínimos quadrados:

> with(linalg):> n:=100;

Warning: new definition for normWarning: new definition for trace

n := 100

Inicialmente usamos uma função para simular os dados experimentais (xi,yi), que sãoarmazenados na matriz A.

> f:=x->exp(x)*sin(2*Pi*x);> A:=[[evalf((2*i-2)/n-1),evalf(f((2*i-2)/n-1))] $i=1..100]:

f := x -> exp(x) sin(2 Pi x)

Aproxima-se agora estes dados com um ajuste dos coeficientes da função base. Testa-seinicialmente a base L={ex, sin x, x, x2, x3}. Uma matriz Bij=Lj(xi), e um vetor vi=yi, sãomontados e então acha-se a solução o através da solução do problema de mínimos quadradosmin [B]{o}-{v}.

> L:=x-> [exp(x),sin(x),x,x^2,x^3];> B:=matrix(n,5,(i,j)->L(A[i][1])[j]):


12

> v:=vector(n,i->A[i][2]): 2 3 L := x -> [exp(x), sin(x), x, x , x ]

> o:=leastsqrs(B,v);

o := [ .1494524053, 1948.933331, -1944.782168, -1.05926618, 305.0439764 ]

> fapp:=x->sum(o[i]*L(x)[i],i=1..5);> print(fapp(x));> plot({A,fapp},-1..1,title=`5 function base`);

5 ----- \ fapp := x -> ) o[i] L(x)[i] / ----- i = 1

2 .1494524053 exp(x) + 1948.933331 sin(x) - 1944.782168 x - 1.05926618 x

3 + 305.0439764 x

Uma vez que a escolha da base foi ineficaz, aumenta-se a ordem dos polinômios atédez., LA={1,x,x2,x3,x4,x5,x6,x7,x8,x9,x10}.

> LA:=x->[x^i $i=0..10];> LA(x);> C:=matrix(n,11,(i,j)->LA(A[i][1])[j]):

i LA := x -> [x $(i = 0 .. 10)]

2 3 4 5 6 7 8 9 10 [1, x, x , x , x , x , x , x , x , x , x ]

> o1:=leastsqrs(C,v);o1 := [ .0015856111, 6.280157300, 6.161396755, -38.14224353, -38.72205101,

60.94847872, 67.27123875, -37.55515556, -46.95928062, 8.475543818,

12.25612434 ]

> fapp2:=x->sum(o1[i]*LA(x)[i],i=1..11);> print(fapp2(x));> plot({A,fapp2},-1..1,title=`10th order polynomial approximation`);


13

11 ----- \ fapp2 := x -> ) o1[i] LA(x)[i] / ----- i = 1

2 3 .0015856111 + 6.280157300 x + 6.161396755 x - 38.14224353 x

4 5 6 7 - 38.72205101 x + 60.94847872 x + 67.27123875 x - 37.55515556 x

8 9 10 - 46.95928062 x + 8.475543818 x + 12.25612434 x

Como (não) se pode ver, a aproximação agora é bem efetiva.

Exemplos de Mecânica dos Sólidos I

Utiliza-se agora o Maple para resolver problemas em Mec. Sol. I. Como primeiroexemplo, constrói-se os diagramas de momento fletor e força cortante da viga em balançoabaixo sujeita ao carregamento mostrado. A força distribuída de 3 kN/m extende-se sobre 8mda viga e uma carga concentrada de 10kN está aplicada sobre a viga curta acoplada.

Usaremos a função pré-programada de Heaviside (degrau unitário) para definir ocarregamento. É fácil notar que o momento fletor é dado por

( )

( )

( ) ( )

M x

Pd xx

x a

Pd a x aa

a x b

Pd a x aa

Pp c Pp x a b b x

( ) =

− ≤

− − +

≤ ≤

− − +

+ − − − ≤

2

2

2

e as forças cortantes por


14

( )( )( )

V x

Pd x x a

Pd a a x b

Pd a Pp b x

( ) =− ≤− ≤ ≤− − ≤

que podem ser escritas em termos de funções de singularidade por> Momento:=(1-Heaviside(x-a))* (-Pd*x*x/2)> +(1-Heaviside(x-b-a))*Heaviside(x-a)* (-Pd*a*(x-a+a/2))> +Heaviside(x-b-a)* (-Pd*a*(x-a+a/2)+Pp*c-Pp*(x-a-b));

Momento := 2 - 1/2 (1 - Heaviside(x - a)) Pd x

- (1 - Heaviside(x - b - a)) Heaviside(x - a) Pd a (x - 1/2 a)

+ Heaviside(x - b - a) (- Pd a (x - 1/2 a) + Pp c - Pp (x - b - a))

> Cortante:=(1-Heaviside(x-a))* (-Pd*x)> +(1-Heaviside(x-b-a))*Heaviside(x-a)* (-Pd*a)> +Heaviside(x-b-a)* (-Pd*a-Pp);

Cortante :=

- (1 - Heaviside(x - a)) Pd x - (1 - Heaviside(x - b - a)) Heaviside(x - a) Pd a

+ Heaviside(x - b - a) (- Pd a - Pp)

Definimos agora um conjunto de dados para traçar um gráfico:> SET:={a=8,b=3,c=2,d=3,Pd=3,Pp=10}:> L:=subs(SET,a+b+c+d);

L := 16

> plot(subs(SET,Momento),x=0..L,title=`BendingMoment`);> plot(subs(SET,Cortante),x=0..L,title=`Shear`);

Reações de Apoio e Deflexão de uma Viga com Carga DistribuídaUma vez que a solução deste problema é simples, pode-se usar diretamente a

integração analítica de equações diferenciais dsolve do Maple. Para resolver-se a reaçãoredundante, utiliza-se o princípio da superposição diretamente, após a solução da viga embalanço.


15

��

��

Po

L

w0

Escreve-se o momento fletor considerando a reação de apoio “Ra” incógnita. Daíresolve-se a equação diferencial da deflexão da viga utilizando-se apenas as condições decontorno da extremidade engastada.d y

dxEI M y L y L

2

20 0= = =( ) , ' ( )

Utiliza-se depois a condição do lado apoiado para descobrir-se a reação desconhecida.> restart;> Momento:=Ra*x-(1/2*(wo*x/L)*x)*x/3;

3 wo x Momento := Ra x - 1/6 ----- L

> deq1:=(D@@2)(y)(x)*EI; (2) deq1 := D (y)(x) EI

> dbound:=y(L)=0,D(y)(L)=0; dbound := y(L) = 0, D(y)(L) = 0

> solution1:=dsolve({deq1=Momento,dbound},y(x));

solution1 :=

5 3 3 2 wo x Ra x L (wo L - 10 Ra) L (wo L - 12 Ra) x y(x) = - 1/120 ----- + 1/6 ----- - 1/30 ----------------- + 1/24 ------------------- EI L EI EI EI

> solve(subs(x=0,rhs(solution1))=0,Ra); 1/10 wo L

Outro problema de deflexão de uma vigaDefine-se o carregamento em termos de funções de Heaviside (degrau unitário) ao

invés de definí-lo como uma função com “if”. Isto evita a integração por intervalo. Obtém-sea solução pela integração das equações diferenciais da viga e apliação das condições decontorno.d y

dxEI M

M dx E I C L

dx y C y L

2

2

1

2

0

0

=

= + =

= + =

∫∫

θ θ

θ

( )

( )


16

��

P

a

L

Define-se a força cortante e obtém-se o momento pela diretamente por integração, uma vezque não há momentos concentrados aplicados. (alternativamente poderia-se definir ocarregamento por funções de Dirac)

> V:=Ra-P*Heaviside(x-L/2); V := Ra - P Heaviside(x - a)

> M:=int(V,x); M := Ra x - (x - a) Heaviside(x - a) P

As rotações são obtidas integrando uma vez o momento. Aplica-se as condições de rotaçãonula no lado engastado para determinar a constante de integração.

> theta:=int(M,x)/EI+C1; 2 2 1/2 Ra x - 1/2 (x - a) Heaviside(x - a) P theta := ------------------------------------------- + C1 EI

> solve(subs(x=L,theta)=0,C1); 2 2 2 Ra L Heaviside(L - a) P L Heaviside(L - a) P L a Heaviside(L - a) P a- 1/2 ----- + 1/2 --------------------- - ---------------------- + 1/2 --------------------- EI EI EI EI

> C1:=%:> print(theta);

2 2 2 2 1/2 Ra x - 1/2 (x - a) Heaviside(x - a) P Ra L Heaviside(L - a) P L ------------------------------------------- - 1/2 ----- + 1/2 --------------------- EI EI EI

2 Heaviside(L - a) P L a Heaviside(L - a) P a - ---------------------- + 1/2 --------------------- EI EI

A seguir integra-se novamente para se obter os deslocamentos. Aplica-se a outra condição decontorno para resolver-se para a outra constante de integração.

> w:=int(theta,x)+C2;

3 3 2 1/6 Ra x - 1/6 (x - a) Heaviside(x - a) P Ra L xw := ------------------------------------------- - 1/2 ------- EI EI

2 2 Heaviside(L - a) P L x Heaviside(L - a) P L a x Heaviside(L - a) P a x + 1/2 ----------------------- - ------------------------ + 1/2 ----------------------- EI EI EI

+ C2

> solve(subs(x=0,w)=0,C2); 3 a Heaviside(- a) P - 1/6 ------------------- EI

> C2:=%:


17

Utiliza-se agora a condição de deslocamento nulo no apoio para obter a reação redundantedesconhecida:

> solve(subs(x=L,w)=0,Ra);

3 2 3 Heaviside(L - a) P L Heaviside(L - a) P L a Heaviside(L - a) P a 3 (1/3 --------------------- - 1/2 ----------------------- + 1/6 --------------------- EI EI EI

3 a Heaviside(- a) P / 3 - 1/6 -------------------) EI / L EI /

> Ra:=%:

Define-se um conjunto de funções para se traçar o gráfico:> SET3:={P=1,L=1,EI=1}:> plot(subs(SET3,w),x=0..1,title=`Deflections`);> plot(subs(SET3,theta),x=0..1,title=`Angle`);> plot(subs(SET3,M),x=0..1,title=`Moment`);


18

Exemplos de Análise Estrutural Avançada

Applicação do Teorema de CastiglianoConsidera-se aqui o mesmo exemplo de uma viga estaticamente indeterminada. Esta

viga está submetica uma carga distribuída triangular e tem um apoio na extremidade esquerdae um engaste na direita. Neste exemplo, obtem-se a reação indeterminada e a rotação no apoioutilizando o segundo teorema de Castigliano. A reação no apoio é obtida da deflexão zero emx=0.

yU

Ra Ra E IM x dx

E IM

M

Radx

M Ra x wx

L

x

L L

o

( ) ( )01

2

10

1

2 3

2

0 0

2

= =

= =

= −

∫ ∫∂∂

∂∂

∂∂

��

�� Po

L

w0

> M:=Ra*x-((1/2)*w0*(x/L)*x)*(x/3); 3 w0 x M := Ra x - 1/6 ----- L

> ya:=int(M*diff(M,Ra)/EI,x=0..L); 3 L (- w0 L + 10 Ra) ya := 1/30 ------------------- EI

> solve(ya=0,Ra); 1/10 w0 L

> a:=Ra=%; a := Ra = 1/10 w0 L

A rotação em x=0 é determinada aplicando-se um momento fictício neste ponto, e depois daintegração este momento é igualado a zero.

> M:=-M0+M; 3 w0 x M := - M0 + Ra x - 1/6 ----- L

> theta0:=subs({M0=0,a},int(M*diff(M,M0)/EI,x=0..L)); 3 L w0 theta0 := - 1/120 ----- EI

Uma Aplicação do Método de Rayleigh-RitzObtem-se a deflexão de uma viga apoiada-engastada, sujeita a uma carga centrada,

caso já resolvido anteriormente por integração direta, mas nesta seção utiliza-se um dos maisimportantes ferramentas da Mecânica, o Método de Rayleigh-Ritz. Calcula-se uma solução


19

aproximada, bastante próxima da solução exata, baseado no princípio da Mínima EnergiaPotencial, que enuncia que a o deslocamento que minimiza a energia potencial de umaestrutura satisfaz as condições de equilíbrio. A Energia Potencial Total consiste na soma daenergia de deformação com o potencial do carregamento aplicado.

( )

F U W

U EId w

dxdx

W P w

L

xL

= +

=

= − −

∫

=

1

2

2

2

2

0

2

A deflexão w no equilíbrio é aquela que minimiza a energia potencial. Aproxima-se a soluçãow por um polinômio. Este polinômio é escolhido de modo a satisfazer as condições decontorno de deflexão nula em ambos os lados e de rotação nula no engaste, tal como

( ) ( ) ( )w x x x L C C x C x= − + +2

1 2 32 , onde os primeiros termos satisfazem as condições de

contorno e o termo final é uma expansão em série de potências truncada no termo quadrático,de modo que há três incógnitas.Busca-se agora os coeficientes C1, C2 e C3 que minimize o potencial F. As condições para a

minimização é a de derivadas primeiras nulas ∂∂

F

C1

0= , ∂∂

F

C2

0= e ∂∂

F

C3

0= .

��

P

a

L

> w:=x*(x-L)^2*(c1+c2*x+c3*x^2); 2 2 w := x (x - L) (c1 + c2 x + c3 x )

> U:=int(EI/2*(diff(w,x,x))^2,x=0..L); 3 2 2 4 2 2 3 U := 2/35 EI L (7 L c2 + 3 L c3 + 35 c1 + 7 L c2 c3)

> W:=P*subs(x=L/2,w); 3 2 W := 1/8 P L (c1 + 1/2 L c2 + 1/4 L c3)

> F:=expand(U+W);F :=

5 2 7 2 3 2 6 2/5 EI L c2 + 6/35 EI L c3 + 2 EI L c1 + 2/5 EI L c2 c3

3 4 5 + 1/8 P L c1 + 1/16 P L c2 + 1/32 P L c3

> solve({diff(F,c1)=0,diff(F,c2)=0,diff(F,c3)=0},{c1,c2,c3}); P P {c2 = - 5/64 ----, c1 = - 1/32 ----, c3 = 0} EI L EI

> w:=subs(%,w); 2 / P P x\ w := x (x - L) | - 1/32 ---- - 5/64 ----| \ EI EI L/

> plot(subs({P=1,L=1,EI=1},w),x=0..1);


20

Segunda Aplicação do Método de Rayleigh-RitzAplica-se o método de Rayleigh-Ritz a uma viga bi-engastada, sujeita a um carregamentouniformemente distribuído sobre metade de seu comprimento. Assim como no caso anterior,utiliza-se uma aproximação polinomical que satisfaz os deslocamentos e rotações nulas em

ambos os lados, neste caso ( ) ( ) ( )w x x x L C C x C x= − + +2 2

1 2 32 . O trabalho das forças

externas é dado neste caso por W Pd wdxL

L= ∫ /2

.

��

L

L/2

��

Pd

��

> w:=x^2*(x-L)^2*(c1+c2*x+c3*x^2+c4*x^3); 2 2 2 3 w := x (x - L) (c1 + c2 x + c3 x + c4 x )

> W:=int(-Pd*w,x=L/2..L); 5 3 2 W := - 1/840 Pd L (28 c1 + 5 c4 L + 8 L c3 + 14 L c2)

5 3 2 + 1/215040 Pd L (3584 c1 + 185 c4 L + 464 L c3 + 1232 L c2)

> U:=int(EI/2*(diff(w,x,x))^2,x=0..L);U := 5 3 4 3 2 1/1155 EI L (165 L c1 c4 + 462 L c1 c2 + 220 L c2 c4 + 297 L c2 c3 + 462 c1

2 2 2 2 6 4 2 5 + 198 L c2 + 264 L c1 c3 + 100 c4 L + 132 L c3 + 220 c4 c3 L )

> F1:=subs({c2=0,c3=0,c4=0},U+W); 5 2 5 F1 := 2/5 EI L c1 - 1/60 Pd L c1

> solve(diff(F1,c1)=0,c1); Pd 1/48 ---- EI

> w1:=subs({c1=%,c2=0,c3=0,c4=0},w); 2 2 x (x - L) Pd w1 := 1/48 -------------- EI


21

> w1:=subs({L=1,EI=1,Pd=1},w1); 2 2 w1 := 1/48 x (x - 1)

> F2:=subs({c3=0,c4=0},U+W); 5 2 2 2 5 F2 := 1/1155 EI L (462 L c1 c2 + 462 c1 + 198 L c2 ) - 1/840 Pd L (28 c1 + 14 L c2)

5 + 1/215040 Pd L (3584 c1 + 1232 L c2)

> solve({diff(F2,c1)=0,diff(F2,c2)=0},{c1,c2}); Pd Pd {c2 = 7/384 ----, c1 = 3/256 ----} EI L EI

> subs(%,subs({c3=0,c4=0},w)); 2 2 / Pd Pd x\ x (x - L) |3/256 ---- + 7/384 ----| \ EI EI L/

> w2:=subs({L=1,EI=1,Pd=1},%); 2 2 w2 := x (x - 1) (3/256 + 7/384 x)

> F3:=U+W;

F3 :=

5 3 4 3 1/1155 EI L (165 L c1 c4 + 462 L c1 c2 + 220 L c2 c4 + 297 L c2 c3

2 2 2 + 462 c1 + 198 L c2

2 2 6 4 2 5 + 264 L c1 c3 + 100 c4 L + 132 L c3 + 220 c4 c3 L )

5 3 2 - 1/840 Pd L (28 c1 + 5 c4 L + 8 L c3 + 14 L c2)

5 3 2 + 1/215040 Pd L (3584 c1 + 185 c4 L + 464 L c3 + 1232 L c2)

>solve({diff(F3,c1)=0,diff(F3,c2)=0,diff(F3,c3)=0,diff(F3,c4)=0},{c1,c2,c3,c4});

Pd 33 Pd 11 Pd 83 Pd {c2 = 1/256 ----, c3 = ---- -----, c4 = - --- -----, c1 = ---- ----} EI L 1024 2 512 3 6144 EI EI L EI L

> subs(%,w); / 2 3\ 2 2 | 83 Pd Pd x 33 Pd x 11 Pd x | x (x - L) |---- ---- + 1/256 ---- + ---- ----- - --- -----| |6144 EI EI L 1024 2 512 3| \ EI L EI L /

> w3:=subs({L=1,EI=1,Pd=1},%); 2 2 / 83 33 2 11 3\ w3 := x (x - 1) |---- + 1/256 x + ---- x - --- x | \6144 1024 512 /

> plot({w1,w2,w3},x=0..1);


22

Exemplo de um Pórtico PlanoAo contrário dos exemplos anteriores, que podem ser resolvidos facilmente pela Mecânicados Sólidos elementar, o pórtico abaixo requer técnicas mais avançadas. Este pórtico consisteem três vigas iguais, com engaste no ponto A e apoio no ponto D, com uma cargauniformemente distribuída P0 aplicada sobre a viga AB. Serão considerados apenas os efeitosde flexão pura e deformação longitudinal; o cisalhamento será desprezado e não há torção.Obter-se-ão as duas componentes da reação de apoio em D utilizando-se o segundo teoremade Castigliano. Observando que não há deslocamentos no ponto F, pode-se escrever que

δ ∂∂x

x

U

R=

*

e δ ∂∂y

y

U

R=

*

, onde a energia de deformação complementar deve ter o mesmo valor

que a energia total de deformação para esta estrutura elástica linear:

UM

EIds

N

EAdsi

i

Li

i

L

i

ni i* = +

∫ ∫∑=

1

2

2

0

2

01

Aqui n é o número total de membros da estrutura, Ni e EAi são as forças e rigidezes, Mi e EIisão os momentos fletores e rigidezes à flexão. Pode-se facilmente notar que as forças axiaisem cada membro são -Ry em CD, Rx em BC e Ry em AB; e os momentos fletores são (-Rx s)em CD, (-Rx L -Ry s) em BC e (-Rx L-Ry L+Rx s+P0 s s/2) em AB.

��

��

L

p0

A

B C

DRx

Ry

> N1:=-Ry;> M1:=-Rx*s;> N2:=Rx;> M2:=-Rx*L-Ry*s;> N3:=Ry;> M3:=-Rx*L-Ry*L+Rx*s+(p0*s)*s/2;


23

N1 := - Ry

M1 := - Rx s

N2 := Rx

M2 := - Rx L - Ry s

N3 := Ry

2 M3 := - Rx L - Ry L + Rx s + 1/2 p0 s

> UN:=(int(N1^2,s=0..L)+int(N2^2,s=0..L)+int(N3^2,s=0..L))/(2*EA);> UM:=(int(M1^2,s=0..L)+int(M2^2,s=0..L)+int(M3^2,s=0..L))/(2*EI);

2 2 2 Ry L + Rx L UN := 1/2 --------------- EA

3 3 3 3 2 3 L (Rx + Ry) Rx L 2 5 3 UM := 1/2 (2/3 Rx L + 1/3 ------------- - 1/3 ------ + 1/20 p0 L + Rx Ry L Ry Ry

4 4 2 3 - 1/12 Rx p0 L - 1/3 p0 Ry L + Ry L )/EI

> deltax:=simplify(diff(UN+UM,Rx)); 2 2 3 L (24 Rx EI + 40 L EA Rx + 24 L EA Ry - L EA p0) deltax := 1/24 --------------------------------------------------- EA EI

> deltay:=simplify(diff(UN+UM,Ry)); 2 2 3 L (12 Ry EI + 6 L EA Rx + 8 L EA Ry - L EA p0) deltay := 1/6 ------------------------------------------------- EA EI

> solve({deltax=0,deltay=0},{Rx,Ry}); 3 2 3 2 L EA p0 (- 3 EI + 4 L EA) L EA p0 (12 EI + 17 L EA) {Rx = - 1/4 --------------------------------, Ry = 1/8 --------------------------------} 4 2 2 2 4 2 2 2 11 L EA + 18 EI + 42 EI L EA 11 L EA + 18 EI + 42 EI L EA

Exemplos em Análise Estrutural Avançada

Faz-se uma análise do erro de interpolação de um elemento finito quadrilateral paraproblemas bidimensional utilizando a expansão em série de Taylor de uma função arbitrária“g” no Maple:

( ) ( ) ( ) ( ) ( )( )g s t s s t ti j

g

s tO s s t t

i

i ji j

i jj

0 00

2

0 0 0

3

0

3

0

2 1,

! != − − + − −

=

+

=∑ ∑ ∂

∂ ∂.

Usa-se o elemento “standard” bilinear isoparamétrico de 4 nós:

( )( ) ( )( )

( )( ) ( )( )

N s t N s t

N s t N s t

1 2

3 4

1

41 1

1

41 1

1

41 1

1

41 1

= − − = + −

= + + = − +

,

,

As primeiras derivadas do erro de interpolação podem serrepresentadas por

( ) ( )∂∂

∂∂

∂∂

∂∂

∂∂

∂∂

e

sg s t

N

s

g

s

e

tg s t

N

t

g

tI

i ii

i

Ii i

i

i

= − = −= =∑ ∑2

1

4

21

4

, , ,

ξ

η

1 2

34

2

2


24

> with(linalg):Warning: new definition for normWarning: new definition for trace

> readlib(mtaylor):

Expansão de Taylor's para a função arbitrária g> gT2:=convert(mtaylor(g(s0,t0),[s0=s,t0=t],3),polynom);

gT2 := g(s, t) + D[1](g)(s, t) (s0 - s) + D[2](g)(s, t) (t0 - t)

2 + 1/2 D[1, 1](g)(s, t) (s0 - s) + (s0 - s) D[1, 2](g)(s, t) (t0 - t)

2 + 1/2 D[2, 2](g)(s, t) (t0 - t)

Define-se um conjunto com os valores em cada nó desta expansão> g1:=subs({s0=-1,t0=-1},gT2):> g2:=subs({s0= 1,t0=-1},gT2):> g3:=subs({s0= 1,t0= 1},gT2):> g4:=subs({s0=-1,t0= 1},gT2):> gv:=[g1,g2,g3,g4];

2 gv := [g(s, t) + D[1](g)(s, t) (- 1 - s) + D[2](g)(s, t) (- 1 - t) + 1/2 %3 (- 1 - s)

2 + (- 1 - s) %2 (- 1 - t) + 1/2 %1 (- 1 - t) ,

2 g(s, t) + D[1](g)(s, t) (1 - s) + D[2](g)(s, t) (- 1 - t) + 1/2 %3 (1 - s)

2 + (1 - s) %2 (- 1 - t) + 1/2 %1 (- 1 - t) ,

2 g(s, t) + D[1](g)(s, t) (1 - s) + D[2](g)(s, t) (1 - t) + 1/2 %3 (1 - s)

2 + (1 - s) %2 (1 - t) + 1/2 %1 (1 - t) ,

2 g(s, t) + D[1](g)(s, t) (- 1 - s) + D[2](g)(s, t) (1 - t) + 1/2 %3 (- 1 - s)

2 + (- 1 - s) %2 (1 - t) + 1/2 %1 (1 - t) ]

%1 := D[2, 2](g)(s, t)

%2 := D[1, 2](g)(s, t)

%3 := D[1, 1](g)(s, t)

Definindo as funções de interpolação e suas derivadas> SF:=[(1-s)*(1-t)/4,(1+s)*(1-t)/4,(1+s)*(1+t)/4,(1-s)*(1+t)/4];> DSFs:=map(x->diff(x,s),SF);> DSFt:=map(x->diff(x,t),SF);

SF :=

[1/4 (1 - s) (1 - t), 1/4 (1 + s) (1 - t), 1/4 (1 + s) (1 + t), 1/4 (1 - s) (1 + t)]

DSFs := [- 1/4 + 1/4 t, 1/4 - 1/4 t, 1/4 + 1/4 t, - 1/4 - 1/4 t]

DSFt := [- 1/4 + 1/4 s, - 1/4 - 1/4 s, 1/4 + 1/4 s, 1/4 - 1/4 s]

A derivada deste erro em relação à coordenada “s” é> deds:=simplify((dotprod(gv,DSFs,'orthogonal'))-D[1](g)(s,t));

deds := - D[1, 1](g)(s, t) s

Similarmente a respeito de “t”> dedt:=simplify((dotprod(gv,DSFt,'orthogonal'))-D[2](g)(s,t));

dedt := - D[2, 2](g)(s, t) t


25

Exemplo de Flexão:Para a geometria descrita pela figura abaixo, considera-se um caso de flexão puracaracterizada pela configuração deformada mostrada com linha sólida na figura. Calcula-se asdeformações deste movimento, assim como as energias de deformação Ix, Iy e Ixynormalizadas em relação a E/(1-ν2). Obtidas estas energias, traça-se um gráfico da razãoIxy/Ix em relação à razão de aspecto hy/hx.

1 2

34

x

y

hx

hy

a

> DSF:=concat(DSFs,DSFt); [ - 1/4 + 1/4 t - 1/4 + 1/4 s ] [ ] [ 1/4 - 1/4 t - 1/4 - 1/4 s ] DSF := [ ] [ 1/4 + 1/4 t 1/4 + 1/4 s ] [ ] [ - 1/4 - 1/4 t 1/4 - 1/4 s ]

> xv:=[0,hx,hx,0]:> yv:=[0,0,hy,hy]:> coor:=concat(xv,yv);

[ 0 0 ] [ ] [ hx 0 ] coor := [ ] [ hx hy ] [ ] [ 0 hy ]

> Jac:=evalm(transpose(coor)*DSF); [ 1/2 hx 0 ] Jac := [ ] [ 0 1/2 hy ]

> Jacinv:=inverse(Jac); [ 2 ] [ ---- 0 ] [ hx ] Jacinv := [ ] [ 2 ] [ 0 ---- ] [ hy ]

> jacobian:=det(Jac); jacobian := 1/4 hx hy

> displace:=matrix(4,2,[-a,0,a,0,-a,0,a,0]); [ - a 0 ] [ ] [ a 0 ] displace := [ ] [ - a 0 ] [ ] [ a 0 ]

> Ddispl:=evalm(transpose(displace)*DSF); [ - a t - a s ] Ddispl := [ ] [ 0 0 ]

> Ddispx:=evalm(Ddispl*Jacinv);


26

[ a t a s ] [ - 2 --- - 2 --- ] Ddispx := [ hx hy ] [ ] [ 0 0 ]

> ex:=Ddispx[1,1];> ey:=Ddispx[2,2];> gxy:=Ddispx[1,2]+Ddispx[2,1];

a t ex := - 2 --- hx

ey := 0

a s gxy := - 2 --- hy

> Ix:=int(int(ex^2*jacobian,s=-1..1),t=-1..1);> Iy:=int(int(ey^2*jacobian,s=-1..1),t=-1..1);> Ixy:=int(int(gxy^2*jacobian,s=-1..1),t=-1..1);

2 a hy Ix := 4/3 ----- hx

Iy := 0

2 a hx Ixy := 4/3 ----- hy

> d33:=subs(nu=3/10,(1-nu)/2); d33 := 7/20

> SET:={hy=x*hx,a=1}:> plot(subs(SET,d33*Ixy/Ix),x=.5..4);

Documents

Uma Introdução ao Maple V - joinville.udesc.br Uma Introdução ao Maple V 2 Álgebra de Matrizes e Vetores A maioria da álgebra linear está incluída no pacote de rotinas do Maple