Exogenous and Endogenous Treatments of Transaction- Transaction-Specific Investments Transaction-Specific Investments

O objeto da aplicação prática será o cálculo da percolação pelas fundações da parede de uma Barragem como mostrado na Figura 1.

Figura 7 – Corte da parede de uma barragem

41

De acordo com a equação (1.16), o fluxo em duas dimensões pela parede de uma barragem de terra pode ser definido por:

� � +

�

� = .

E no caso de um meio isotrópico, temos a Equação de Laplace em duas dimensões:

� � +

�

� = .

A modelagem da barragem apresentará uma simplificação para facilitar o cálculo da Solução Analítica, tendo em vista que o objetivo do estudo é a comparação entre esta e a Solução Numérica obtida através do Método sem malhas.

A simplificação consiste em se considerar apenas o fluxo que acontece abaixo da parede da barragem, representado pelas linhas na figura. O obstáculo é o material que forma o núcleo impermeável dessa barragem de terra, portanto, o fluxo acontecerá apenas no maciço dessa barragem. Ressaltando novamente que este é um caso didático onde se objetiva apenas a aplicação da rotina em questão, mas que na prática da engenharia não é válido.

Fica estabelecido o seguinte domínio como mostrado na Figura 7:

 15 metros abaixo do maciço da parede – eixo x.

 20 metros ao longo da seção transversal dessa parede – eixo y.

Como o fluxo é nulo em u(x,0), u(0,y), u(20,y), temos as seguintes condições de contorno:

, = .

, = .

, = .

42

Para se encontrar a condição de contorno em u(x,15), calcula-se inicialmente o gradiente hidráulico de percolação de acordo com (1.2):

=Δu= − = , / .

Dessa forma temos como condição de contorno em u(x,15):

, = − , .

Temos agora as seguintes condições de fronteira:

, = .

, = .

, = .

, = − , .

A partir do valor dessas condições de fronteira, calcula-se através de Série de Fourier o coeficiente cn através da equação (1.56):

=

ℎ � .

Aplicando esse coeficiente na equação (1.55) temos então a Solução Analítica dada por:

, =

� _{ℎ �}

43

2.8 Processamento da Rotina

O Scilab tem comportamento semelhante ao Matlab e outros pacotes computacionais de grande capacidade de processamento. Foi desenvolvido por pesquisadores franceses desde o início da década de 90, ele é distribuído de forma gratuita desde 1994, e por esse motivo vem tendo seu uso bastante difundido em pesquisas e no meio acadêmico.

A Rotina utilizada se encontra no Anexo A. Na primeira tela são solicitados os dados da equação diferencial:

Figura 8 – Valores dos coeficientes da EDP.

Fonte: elaborado pelo autor

44

Figura 9 _{– Número de pontos, e dimensões do domínio.}

Fonte: elaborado pelo autor

Em seguida vem a janela onde deverão ser inseridas as condições de contorno do problema:

Figura 10 – Condições de contorno.

Fonte: elaborado pelo autor

Em seguida é dada a opção de se informar o parâmetro de forma C, ou fazer a busca por um C ótimo:

45

Figura 11 _{– Escolha do parâmetro de forma C.}

Fonte: elaborado pelo autor

Como iremos buscar a otimização, a seguir pergunta-se qual a precisão, ou opção de refinamento desse parâmetro de forma:

Figura 12 – Refinamento do parâmetro de forma.

Fonte: elaborado pelo autor

A partir desse ponto o programa roda a rotina, fazendo a busca pelo parâmetro de forma ótimo.

46

2.9 Resultados

Foram realizadas as seguintes simulações:

Com um refinamento de 10-1_{, foram realizadas simulações com:}

 Discretização com 9 pontos  Discretização com 16 pontos  Discretização com 25 pontos

Com um refinamento de 10-2_{, foram realizadas simulações com:}

 Discretização com 9 pontos  Discretização com 16 pontos  Discretização com 25 pontos

Com um refinamento de 10-3_{, foram realizadas simulações com:}

 Discretização com 9 pontos  Discretização com 16 pontos

Os gráficos inicialmente mostram a discretização dos pontos no domínio. Os pontos foram distribuídos linearmente através do domínio. Abaixo estão os gráficos dos pontos ao longo do domínio:

47

Figura 13 _{– Discretização com 9 pontos.}

Fonte: elaborado pelo autor

Figura 14 – Discretização do domínio com 16 pontos.

48

Figura 15 – Discretização do domínio com 25 pontos.

Fonte: elaborado pelo autor

Para se conseguir um universo mais amplo de situações, foi adotado um intervalo para o valor de C de 0,01 a 1000. Devido a esse intervalo ser bem grande, utilizou-se uma escala logarítmica no gráfico para que ficasse demonstrado o comportamento do parâmetro em todo o intervalo amostral.

Em cada gráfico do fluxo, estão a Solução Numérica e Analítica para uma melhor visualização.

Foram realizadas simulações utilizando a busca pelo parâmetro C otimizado, e simulações e com o C calculado através da fórmula (2.6):

= √ �

√ ��−

49

Inicialmente foi realizada a simulação com os valores calculados pela equação (2.6), de onde obtivemos para cada situação os seguintes valores de C:

= √ � √ ��− =√ × √ − = √ − = , . = √ � √ ��− =√ × √ − = √ − = , . = √ � √ ��− =√ × √ − = √ − = , .

Onde AS é a área total do domínio, e _�� é o número de pontos do domínio.

50

Figura 16 – Carga hidráulica u(m) considerando C estimado (C=8,66) e 9 pontos no domínio.

51

Figura 17 – Carga hidráulica u(m) considerando C estimado (C=5,77) e 16 pontos no domínio.

52

Figura 18 – Carga hidráulica u(m) considerando C estimado (C=4,33) e 25 pontos no domínio.

53

Figura 19 – Carga hidráulica u(m) considerando C otimizado (C=449,2) e 9 pontos no domínio, com refinamento 10

-1

_.

54

Figura 20 – Resíduos com C otimizado: 9 pontos com refinamento 10

-1

_.

55

Figura 21 – Carga hidráulica u(m) considerando C otimizado (C=436,7) e 16 pontos no domínio, com refinamento 10

-1

_.

56

Figura 22 – Resíduos com C otimizado: 16 pontos com refinamento 10

-1

_.

57

Figura 23 – Carga hidráulica u(m) considerando C otimizado (C=79,2) e 25 pontos no domínio, com refinamento 10

-1

_.

58

Figura 24 – Resíduos com C otimizado: 25 pontos com refinamento 10

-1

_.

59

Figura 25 – Carga hidráulica u(m) considerando C otimizado (C=0,01) e 9 pontos no domínio, com refinamento 10

-2

_.

60

Figura 26 – Resíduos com C otimizado: 9 pontos com refinamento 10

-2

_.

61

Figura 27 – Carga hidráulica u(m) considerando C otimizado (C=17,46) e 16 pontos no domínio, com refinamento 10

-2

_.

62

Figura 28 – Resíduos com C otimizado: 16 pontos com refinamento 10

-2

_.

63

Figura 29 – Carga hidráulica u(m) considerando C otimizado (C=80,23) e 25 pontos no domínio, com refinamento 10

-2

_.

64

Figura 30 – Resíduos com C otimizado: 25 pontos com refinamento 10

-2

_.

65

Figura 31 – Carga hidráulica u(m) considerando C otimizado (C=0,001) e 9 pontos no domínio, com refinamento 10

-3

_.

66

Figura 32 – Resíduos com C otimizado: 9 pontos com refinamento 10

-3

_.

67

Figura 33 – Carga hidráulica u(m) considerando C otimizado (C=17,464) e 16 pontos no domínio, com refinamento 10

-3

_.

68

Figura 34 – Resíduos com C otimizado: 16 pontos com refinamento 10

-3

_.

69

2.10 Análise de Resultados

As simulações realizadas encontraram diferentes valores de C para cada situação, onde variaram o número de pontos discretizados domínio, e a precisão do parâmetro de forma.

Abaixo segue um quadro comparativo dos valores encontrados de C:

Tabela 13 – Comparativo dos valores encontrados de C.

Pontos discretizados Valor de C estimado Valores otimizados de C = √�� √���− Refinamento de 10-1 Refinamento de ₁₀-2 Refinamento de ₁₀-3 9 8,66 449,2 0,01 0,001 16 5,77 436,7 17,46 17,464 25 4,33 79,2 80,23

Fonte: elaborado pelo autor

Para o refinamento de 10-1_{, o intervalo de busca utilizado para se encontrar o parâmetro C foi de 0 a}

1000. Ao se aplicar esse mesmo intervalo para o refinamento de 10-2_{, o tempo de processamento dos}

dados foi muito extenso, sendo que para o caso de apenas 9 pontos chegou a quase duas horas. Dessa forma optou-se por fazer as demais simulações no intervalo de 0 a 100, o que também apresentou resultados condizentes com os observados nos gráficos com intervalo maior.

Observou-se uma discrepância entre os valores encontrados para o parâmetro de forma calculado através da fórmula (2.6), e os valores encontrados através da otimização. Ao observar-se os gráficos nota-se que as simulações em que se utilizou o C estimado pela equação (2.6), apresentou uma aproximação inferior às obtidas com os gráficos obtidos com o C otimizado pela rotina.

70

Nas simulações tivemos 2 características variáveis:

 Número de pontos discretizados no domínio  Refinamento da busca

Fazendo uma análise de cada uma dessas características observa-se que:

 Quanto mais pontos utilizar-se na discretização do domínio, maior a aproximação dos gráficos da solução analítica e solução numérica.

 O refinamento no valor do parâmetro não se mostrou tão importante quanto o número de pontos, pois analisando o gráfico da Figura 23 (25 pontos, refinamento 10-1_{), e o}

gráfico da Figura 31 (9 pontos, refinamento 10-3_{), observa-se que o primeiro apresenta}

mais pontos, e uma maior aproximação entre a Solução Numérica e a Solução Analítica, enquanto o segundo apresenta um melhor refinamento, mas não apresenta uma aproximação tão boa entre a Solução Numérica e a Solução Analítica.

71

3. Conclusões

O que buscou-se nesse estudo foi desenvolver uma rotina que calculasse o parâmetro C otimizado para um determinado intervalo, e utilizar esse parâmetro no método numérico sem malhas, Meshless, para obter uma solução numérica para uma equação diferencial parcial.

A rotina desenvolvida forneceu as simulações apresentadas anteriormente, e o que se observou nessas simulações é que a precisão da Solução Numérica aumenta com a quantidade de pontos, ou seja, quanto maior o número de pontos discretizados no domínio, melhor será a aproximação entre a Solução Numérica e a Solução Analítica. Em relação ao refinamento não ocorreu o mesmo, já que nos resultados analisados, por mais que se aumente a quantidade de casas decimais do C, não teremos um aumento tão significativo na aproximação entre a Solução Numérica e Solução Analítica.

Embora o desempenho da rotina possa ser considerado bom para o caso prático proposto, vela ressaltar que em muitas aplicações da Equação de Laplace, seria interessante se usar um refinamento de no mínimo 10-5 _{(Colier et al., 1997), o que não foi possível no nosso caso, onde chegamos até a 10}-3_.

Portanto, conclui-se que no caso estudado, em que aplicamos uma equação diferencial parcial para determinar o fluxo através da parede de uma barragem, o resultado poder ser considerado satisfatório, tendo em vista que como dito no capítulo sobre barragens, é muito importante se determinar o caminho seguido pelo fluxo, ou seja, o percurso que a água irá percorrer dentro do maciço da barragem nesse processo de percolação, e o que notou-se foi uma boa aproximação geométrica, visto que nos casos em que discretizou-se o domínio com 25 pontos, a Solução Analítica e a Solução Numérica são bem próximas.

Como uma sugestão para um trabalhos futuros seguem algumas indicações:

 Complementação da rotina, de modo a calcular o erro absoluto para as simulações obtidas;

72

 Aplicação de outros tipos de Funções de Base Radial (nesse trabalho só foi explorado

a RBF multiquadrática);

73

REFERÊNCIAS

Belinha, J. (2010). The Natural Neighbour Radial Point Interpolation. 2010 Dissertação (Doutoramento em Engenharia Mecânica) - Faculdade de Engenharia da Universidade do Porto,Porto 2010.

Belytschko, T. ; Liu Y.Y., Gu. L. (1994) Element Free-Galerkin. International for Numerical Methods in Fluids, vol. 37: 229-256.

Borges, Altemir José (2006), “Notas de Aula da disciplina Equações Diferenciais Parciais”. Curso de Bacharelado em matemática da UFTPR, Curitiba.

BOYCE, W. E. e PRIMA, Richard C. Di. Equações diferenciais elementares e problemas de valores de contorno. RJ, LTC, 1994.

BRONSON, Richard. Moderna introdução às equações diferenciais. SP: McGraw-Hill, 1977.

BUHMANN, M. D., (2003), Radial Basis Function: Theory and Implementations, Cambridge University Press, Cambridge.

CAPUTO, HOMERO PINTO. Mecânica dos solos e suas aplicações, volume 1: fundamentos. 6 ed. Rio de Janeiro: LTC, 2007.

CHENG, A. H. D., GOLBERG, M. A., KANSA, E. J., ZAMMITO, G., (2003), “Exponential convergence and H-c multiquadric collocation method for partial differential equations”, Numerical Methods for Partial Differential Equations, 19(5) : 571-594.

Cingoski, V.; Miyamoto, N. ; Kaneda, K. and Yamashita, H. (1998). Element-Free Galerkin Method for Electromagnetic Field Computations. IEEE Transactions on Magnetics, vol. 34, n. 5, September 1998.

Colier, B.; Simkin, J.; Trowbridge, C. W.; Barberis, U.; Picco, E.; Gutierrez, T.; Longo, A.; Greenough, C.; Thomas, D.; Alloto, P.; Molfino, P.; Molinari, G. (1997) Project MIDAS: Magnet Integrated Design and Analysis System, IEEE Transactions on Magnetics, 33, no 2, March, 1143- 1148.

Cruz, P. T. (1996) “Permeabilidade e condutividade”, in 100 barragens brasileiras: casos históricos, materiais de construção, projeto. Oficina de Textos, São Paulo – SP, pp. 258 – 278.

Forger, Michael (2012), “Notas de aula da disciplina Equações Diferenciais Parciais”. Curso de Mestrado em Matemática da USP, São Paulo.

HON, Y. C., CHEUNG, K. F., MAO, X. Z., KANSA, E.J., (1999), “Multiquadric solution for shallow water equations”, Journal of Hydraulic Engineering-ASCE, 125(5) : 524-533.

74

ICOLD Bulletin 128, Management of Reservoir Water Quality - Introduction and Recommendations, 2004.

ICOLD, Technical Dictionary on Dams - Glossary of Terms, 1994.

KANSA, E.J., (1990), “Multiquadrics - a scattered data approximation scheme with applications to computational fluid dynamics II: solutions to parabolic, hyperbolic and elliptic partial differential equations”, Computers and Mathematics with Application.

KANSA, E. J., (1992), “A strictly conservative spatial approximation scheme for the governing engineering and physics equations over irregular regions and inhomogeneaus scattered nodes”, Computers and Mathematics with Applications, 24(5-6) : 169-190.

KANSA, E. J., HON, Y.C., (2000), “Circumventing the ill-conditioning problem with multiquadric radial basis functions: applications to elliptic partial differential equations”, Computers and Mathematics with Applications, 39(7-8) : 123-137.

León, L. A. M. ; CASTRO, R. G. S. (2010) Equações Diferenciais Parciais. Notas de aula do Curso de Matemática da Universidade Estadual do Norte Fluminense, Rio de Janeiro.

Liu, Gui-Rong. (2003), Mesh free methods : moving beyond the finite element method, CRC, New York.

Liu, W. K. ; Jun, S. and Zhang, F. (1995 a) Reproducing Kernel Particle Method. International for Numerical Methods in Fluids, vol. 20: 181-1106.

Liu, W.K.; S. Li, and T. Belytschko. (1995b) Moving Least Square Reproducing Methods Methodology and Convergence. Technical Report 95-3-XX, TICAM, The University of Texas at Austin, 1995.http://mecad.uta.edu/~jinman/meshless.html.

Monaghan, J.J. (1982) Smooth Particle Hydrodynamics. Anual Review of Astronomy and Astrophysics. vol.30: 543-574.

Nayroles, B. G. Touzot, and P. Villon. (1992) Generalizing the Finite Element Method: Diffuse Approximation and Diffuse Element. Computational Mechanics, vol. 10: 337-318.

SANDRONI, S., (2012), “Notas de aula da disciplina de Barragens de Terra e Enrocamento”. Curso de Mestrado da COPPE/UFRJ, Rio de Janeiro.

Santos, Reginaldo J. (2007), “Notas de aula da disciplina Equações Diferenciais Parciais”. Curso de bacharelado em Matemática da UFMG, Rio de Janeiro.

Santos, Reginaldo J. (2011), Introdução às Equações Diferenciais Ordinárias, Imprensa Universitária da UFMG, Belo Horizonte, 2011.

SCHABACK, R., WENDLAND, H., (2002), “Approximation by positive definite kernels”, IDOMat2001 proceedings (invited lecture).

75

SOUZA, V. A. D. (2005) Simulação do regime de Fluxo no Maciço de Terra Compactada da Barragem Jaburu I – Dissertação de Mestrado, Universidade Federal do Ceará, Fortaleza – CE, pp. 33 – 34.

Viana S. A. ; Mesquita, R.C. (1998a) Moving Least Square Reproducing Kernel Method for Electromagnetic Field Computation. The Eighth Biennial IEEE Conference on Electromagnetic Field Computation, Tucson, Arizona, Junho de 1998.

WENDLAND, H., (1995), “Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree”, Advances in Computational Mathematics, 4 : 389-396.

76

ANEXO A – ROTINA EM SCILAB

//========================================================================= =====

//...DISCIPLINA: MÉTODOS NUMÉRICOS

//...CURSO: MESTRADO EM RECURSOS HÍDRICOS //...UNIVERSIDADE FEDERAL DO CEARÁ (UFC)

//...CÁLCULO DA PERCOLAÇÃO PELAS FUNDAÇÕES DE UMA BARRAGEM USANDO MESHLESS //...PROF.: MARCO AURÉLIO HOLANDA DE CASTRO

//...ALUNO HÉRCULES LIMA DE MEDEIROS //...AGOSTO DE 2014

//========================================================================= =====

clear(); //Apaga todas as variáveis da memória

txt=['A=';'B=';'C=';'D=';'E=';'F=';'G=';'H=';'I='];//Constantes EDO

valor=x_mdialog('Eq. Geral: A.∂²u/∂x² + B.∂²u/∂y² + C.∂u/∂x + D.∂u/∂y + E.u(x) + F.u(y) + G.x + H.y = I.f(x,y)',txt,['1';'1';'0';'0';'0';'0';'0';'0';'0']); //Sugere valores

A=evstr(valor(1)); B=evstr(valor(2)); C=evstr(valor(3)); D=evstr(valor(4)); E=evstr(valor(5)); F=evstr(valor(6)); G=evstr(valor(7)); H=evstr(valor(8)); I=evstr(valor(9));

//========OBTENDO VALORES DO PROBLEMA================= txt1=['Número de Pontos';'Lx=';'Ly='];//Constantes EDO

valor1=x_mdialog('Lados do Retângulo (m):',txt1,['25';'20';'15']); //Sugere valores

NP=evstr(valor1(1)); Lx=evstr(valor1(2)); Ly=evstr(valor1(3)); rNP=int(NP^0.5) dx=Lx/(rNP-1) dy=Ly/(rNP-1) x=0 k=0 for i=1:(rNP) if i<>1 then if i<>(rNP) then x=x+dx y=0 for j=1:(rNP) if j<>1 then if j<>(rNP) then k=k+1 y=y+dy Xi(k)=x Yi(k)=y fxy(k)=I end end end end end end x=0 [o l]=size(Xi)

77

for i=1:(rNP) y=0 for j=1:(rNP) if i==1 then if j==1 then k=k+1 Xi(k)=0 Yi(k)=0 else k=k+1 y=y+dy Xi(k)=0 Yi(k)=y end else if j==1 then k=k+1 Xi(k)=x Yi(k)=0 else y=y+dy if i==rNP then k=k+1 Xi(k)=x Yi(k)=y else if j==rNP then k=k+1 Xi(k)=x Yi(k)=y end end end end end x=x+dx end x=Xi y=Yi

txt2=['Condições de Contorno: f(0,y)=';'f(x,0)=';'f(Lx,y)=';'f(x,Ly)='];//Constantes EDO valor2=x_mdialog('Equações do contorno:',txt2,['0';'0';'0';'20-0.25*x']); //Sugere valores

F0y=evstr(valor2(1));

Fx0=evstr(valor2(2));

Fny=evstr(valor2(3));

Fxn=evstr(valor2(4)); k=o for i=1:(rNP) for j=1:(rNP) if i==1 then if j==1 then k=k+1

if sum(F0y)<>0 then fxy(k)=F0y(k) else fxy(k)=0 end else k=k+1

if sum(F0y)<>0 then fxy(k)=F0y(k)

78

else fxy(k)=0 end end else if j==1 then k=k+1 if sum(Fx0)<>0 then fxy(k)=Fx0(k) else fxy(k)=0 end else if i==rNP then k=k+1

if sum(Fny)<>0 then fxy(k)=Fny(k) else fxy(k)=0 end else if j==rNP then if sum(Fxn)<>0 then k=k+1

fxy(k)=400*(sinh(3.141592654*y(k)/20)*sin(3.141592654*x(k)/20))/sinh(3.141592654*15/20)

else fxy(k)=0 end end end end end end end Xj=Xi Yj=Yi plo=[x y fxy] disp(plo)

tic()//computa o tempo de cpu

//===========================CALCULO DA MATRIZ R================================

[m n]=size(Xi)

for i=1:m //Laço para linhas for j=1:m //Laço para as colunas

R(i,j)=sqrt([Xi(i,1)-Xj(j,1)]^2+[Yi(i,1)-Yj(j,1)]^2); //Determina a matriz R end

end

for i=1:m //Laço para as Linhas

Man(i,1)=400*(sinh(3.141592654*Yi(i)/20)*sin(3.141592654*Xi(i)/20))/sinh(3.141592654*15/20)

end

//=============OBTEM INFORMAÇÕES DO PARÂMETRO FORMA=============================

[num]=x_choose(['Informar o parâmetro de forma';'Obter o parâmetro de forma ótimo'],...

['Selecionar a opção com duplo clique'],'Cancela')

if num==2|num==0 then

labels=["Refinamento do Parâmetro C:"];

79

list("vec",1),["10^-1"])

//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +++++

//---INICIA A BUSCA DO VALOR ÓTIMO DE C DA EDO ---

//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +++++

//=====================INICIALIZA

VARIÁVEIS===================================== [m n]=size(Xi);// Número de pontos do problema

PF=vref; //Parâmetro Inicial a ser testado segue=0; // varável de repetição

p=0 src=0 srd=0 sri=0//log(ref) sra=0 conta=0 inicia=0 mpf1=0 srda1=0 srda2=0 k=0 while PF<=1000 k=k+1 sra=sri srca=src srda=srd //==========================MONTA A MATRIZ FI=================================== [m n]=size(R)//obtem a quantidade de linha (ri) e colunas (rj) for i=1:m //Laço para linhas

for j=1:n //Laço para as colunas

FI(i,j)=(R(i,j)^2+PF^2)^0.5; //Determina a matriz FI end //end for j

end// end for i

//==================CALCULA A MATRIZ

L(FI)====================================== [m n]=size(FI); //descobre a dimensão da matriz FI

//Determina L(FI) para 1ª derivada em X for i=1:m //Laço para linhas

for j=1:n //Laço para as colunas

LFIx(i,j)=(Xi(i)-Xj(j))/sqrt(R(i,j)^2+PF^2)

end//end for j end //end for i

//Determina L(FI) para a 1ª derivada em Y for i=1:m //Laço para linhas

for j=1:n //Laço para as colunas

LFIy(i,j)=(Yi(i)-Yj(j))/sqrt(R(i,j)^2+PF^2)

end//end for j end //end for i

//Determina L(FI) para a 2ª derivada em X for i=1:m //Laço para linhas

for j=1:n //Laço para as colunas

LFIx2(i,j)=((Yi(i)-Yj(j))^2+PF^2)/(R(i,j)^2+PF^2)^1.5

end//end for j end //end for i

//Determina L(FI) para a 2ª derivada em Y for i=1:m //Laço para linhas

80

for j=1:n //Laço para as colunas

LFIy2(i,j)=((Xi(i)-Xj(j))^2+PF^2)/(R(i,j)^2+PF^2)^1.5

end//end for j end //end for i

LFI=(A*LFIx2)+(B*LFIy2)+(C*LFIx)+(D*LFIy); //Monta a matriz LFI de acordo com a equação //===================DETERMINA A MATRIZ COLOCAÇÃO

'A'============================ A1=LFI(1:o,1:m); //recorta Linha 1 até Npontos A2=FI((o+1):m,1:m);//Recorta a linha de nPontos até m MA=[A1;A2]; //Concatena a matriz A1 com A2

VetorY=fxy //vetorY sem considerar as condições de contorno

//==================DETERMINA A MATRIZ LAMBIDA - ML============================ ML=MA\VetorY;

//=========CALCULA O RESÍDUO PARA DETERMINAR O "C" ÓTIMO (PF)================== res=0;

GXY=FI*ML;

srd=log10(sum(abs(abs(GXY(1:o))-abs(Man(1:o))))/o);

src=log10(sum(abs(abs(GXY(o+1:m))-abs(Man(o+1:m))))/(m-o)); plo=[src srd, PF] disp(plo) mpf=PF srda1=srd srca1=src Srd(k)=srd Src(k)=src Mpf(k)=PF if srda2==0 then srda2=srda1 srca2=srca1 mpf1=PF end

if srda1<srda2 then srda2=srda1 srca2=srca1 mpf1=PF end PF=PF+vref end if mpf1 <>0 then PF=mpf1 else PF=mpf end end //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +++++

//---TERMINA A OTIMIZAÇÃO DE 'C' AQUI---

//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +++++

if num==1 then

labels=["Parâmetro de Forma (C):"];

[ok,PF]=getvalue("Informar o parâmetro de forma(C)",labels,...

list("vec",1),["7.01"])

//PF=PFindicado end

//---MONTA O VETOR Y---

81

//MONTA A MATRIZ FI

[m n]=size(R)//obtem a quantidade de linha (ri) e colunas (rj) for i=1:m //Laço para linhas

for j=1:n //Laço para as colunas

FI(i,j)=(R(i,j)^2+PF^2)^0.5; //Determina a matriz FI end //end for j

end// end for i

//---CALCULA A MATRIZ L(FI)---

[m n]=size(FI); //descobre a dimensão da matriz FI //Determina L(FI) para 1ª derivada em X

for i=1:m //Laço para linhas for j=1:n //Laço para as colunas

LFIx(i,j)=(Xi(i)-Xj(j))/sqrt(R(i,j)^2+PF^2)

end//end for j end //end for i

//---Determina L(FI) para a 1ª derivada em Y--- for i=1:m //Laço para linhas

for j=1:n //Laço para as colunas

LFIy(i,j)=Yi(i)-Yj(j)/(sqrt(R(i,j)^2+PF^2))

end//end for j end //end for i

//---Determina L(FI) para a 2ª derivada em X--- for i=1:m //Laço para linhas

for j=1:n //Laço para as colunas

LFIx2(i,j)=((Yi(i)-Yj(j))^2+PF^2)/(R(i,j)^2+PF^2)^1.5

end//end for j end //end for i

//---Determina L(FI) para a 2ª derivada em Y--- for i=1:m //Laço para linhas

for j=1:n //Laço para as colunas

LFIy2(i,j)=((Xi(i)-Xj(j))^2+PF^2)/(R(i,j)^2+PF^2)^1.5

end//end for j end //end for i

LFI=(A*LFIx2)+(B*LFIy2)+(C*LFIx)+(D*LFIy); //Monta a matriz LFI de acordo com a equação //---DETERMINA A MATRIZ COLOCAÇÃO 'A'---

A1=LFI(1:o,1:m); //recorta Linha 1 até Npontos A2=FI((o+1):m,1:m);//Recorta a linha de nPontos até m MA=[A1;A2]; //Concatena a matriz A1 com A2

//---DETERMINA A MATRIZ LAMBIDA - ML--- ML=MA\VetorY

//ML=MA\VetorY; //produto da matriz A inversa pelo vetor y

//---DETERMINA A MATRIZ DE APROXIMAÇÃO - Map = FI*ML--- Map=FI*ML; //matriz com os pontos de aproximação (SOLUÇÃO)

//---Determina matriz solução analítica--- M=[Xi Yi Map]

if getos()=='Windows' then

unix('del scotm_cc.txt');//arquivo saída else unix('rm -f scotm_cc.txt'); end write('scotm_cc.txt',M); // ========================================================================== ===

82

//Usa um script de interpolação do toolbox surface (www.reveyrand.fr) para... //espacializar o resultado no domínio. Autor: Tibault Reveyrand

//

========================================================================== ===

functionoutput=interpolate_grid(M, scaling_factor)//interpolate_grid X_min=min(M(:,1)); X_max=max(M(:,1)); Y_min=min(M(:,2)); Y_max=max(M(:,2)); Z_min=min(M(:,3)); Z_max=max(M(:,3));

step=min(abs((M(1:$-1,1)-M(2:$,1))+%i*(M(1:$-1,2)-M(2:$,2))))/scaling_factor; X_nbp=1+int((X_max-X_min)/step);

Y_nbp=1+int((Y_max-Y_min)/step); x=linspace(0,1,X_nbp).';//X_nbp y=linspace(0,1,Y_nbp).';//Y_nbp z=zeros(X_nbp,Y_nbp);

M(:,1)=(M(:,1)-X_min)./(X_max-X_min); M(:,2)=(M(:,2)-Y_min)./(Y_max-Y_min); M(:,3)=(M(:,3)-Z_min)./(Z_max-Z_min); for i=1:X_nbp,

for j=1:Y_nbp,

d=abs(M(:,1)-x(i)+%i*(M(:,2)-y(j)))./max(abs(M(:,1)-x(i)+%i*(M(:,2)-y(j))));

w=(((1-(d)).^2).*(1+2.*(d)))./(1-(((1-(d)).^2).*(1+2.*(d)))+((1-(((1-(d)).^2).*(1+2.*(d))))==0)*%eps); MT=M(find(w~=0),:);

w=w(find(w~=0));

_system=[w.*MT(:,1).^2,w.*MT(:,2).^2,w.*MT(:,1).*MT(:,2),w.*MT(:,1),w.*MT(:,2),w.*((MT(:,1).*0)+1)]; _coeff=pinv(_system)*(w.*MT(:,3));

z(i,j)=(_coeff(1)*x(i)^2+_coeff(2)*y(j)^2+_coeff(3)*x(i)*y(j)+_coeff(4)*x(i)+_coeff(5)*y(j)+_coeff(6)); end;

end;

output=list(x.*(X_max-X_min)+X_min,y.*(Y_max-Y_min)+Y_min,z.*(Z_max-Z_min)+Z_min);

endfunction

z=interpolate_grid(M,10); N=[Xi Yi Man]

if getos()=='Windows' then

unix('del scotm_cc.txt');//arquivo saída else unix('rm -f scotm_cc.txt'); end //write('scotm_cc.txt',M); // ========================================================================== ===

//Usa um script de interpolação do toolbox surface (www.reveyrand.fr) para... //espacializar o resultado no domínio. Autor: Tibault Reveyrand

//

========================================================================== ===

functionoutput=interpolate_grid(N, scaling_factor)

X_min=min(N(:,1)); X_max=max(N(:,1)); Y_min=min(N(:,2)); Y_max=max(N(:,2)); Z_min=min(N(:,3)); Z_max=max(N(:,3));

step=min(abs((N(1:$-1,1)-N(2:$,1))+%i*(N(1:$-1,2)-N(2:$,2))))/scaling_factor; X_nbp=1+int((X_max-X_min)/step);

Y_nbp=1+int((Y_max-Y_min)/step); x=linspace(0,1,X_nbp).';

83

u=zeros(X_nbp,Y_nbp);

N(:,1)=(N(:,1)-X_min)./(X_max-X_min); N(:,2)=(N(:,2)-Y_min)./(Y_max-Y_min); N(:,3)=(N(:,3)-Z_min)./(Z_max-Z_min); for i=1:X_nbp,

for j=1:Y_nbp,

d=abs(N(:,1)-x(i)+%i*(N(:,2)-y(j)))./max(abs(N(:,1)-x(i)+%i*(N(:,2)-y(j))));

w=(((1-(d)).^2).*(1+2.*(d)))./(1-(((1-(d)).^2).*(1+2.*(d)))+((1-(((1-(d)).^2).*(1+2.*(d))))==0)*%eps); NT=N(find(w~=0),:);

w=w(find(w~=0));

_system=[w.*NT(:,1).^2,w.*NT(:,2).^2,w.*NT(:,1).*NT(:,2),w.*NT(:,1),w.*NT(:,2),w.*((NT(:,1).*0)+1)]; _coeff=pinv(_system)*(w.*NT(:,3));

u(i,j)=(_coeff(1)*x(i)^2+_coeff(2)*y(j)^2+_coeff(3)*x(i)*y(j)+_coeff(4)*x(i)+_coeff(5)*y(j)+_coeff(6)); end;

end;

output=list(x.*(X_max-X_min)+X_min,y.*(Y_max-Y_min)+Y_min,u.*(Z_max-Z_min)+Z_min);

endfunction

u=interpolate_grid(N,10); //Atribui a z os valores interpolados

//Obs: caso queira modificar o espaçamento da grade, editar o fator de escala (M,fator)

//===============APRESENTA O TEMPO DE PROCESSAMENTO============================= tempo_exec=toc()

disp('Tempo de execução (s)= '+string(tempo_exec))

//==========IMPRESSÃO DE

RESULTADOS============================================= if getos()=='Windows' then

unix('del sz3cotm_cc.txt');//arquivo saída else

unix('rm -f sz3cotm_cc.txt'); end

//teste=[Map Man Man-Map]

disp(Map) disp(Man) disp(Man-Map)

write('sz3cotm_cc.txt',u(3));

scf()//Cria nova janela para impressão do gráfico plot(Xi,Yi,'o'); //Plota os pontos do domínio

In document Creating the Image: A Transaction Cost Analysis of Joint Value Creation in the Motion Picture Industry (sider 52-58)