EXERCÍCIO COMPUTACIONAL – MAP 5725 RELATÓRIO · mÉtodo de runge-kutta-fehlberg (rkf) O método RKF utiliza os métodos RK45 e RK56, a idéia desse método é obter uma estimativa

1

EXERCÍCIO COMPUTACIONAL – MAP 5725

RELATÓRIO

NOME: Rafael Fernandes Pinheiro

N.USP.: 5383783

QUESTÃO 1

Queremos resolver: �� (I)

Para Resolver essa equação usamos a seguinte transformação: �� (II)

derivando em relação a �, temos: ��

�� (III)

Substituindo (I) em (II), temos: ��

�� 1 � �� .

Obtemos: �� 1 � �� (IV).

Essa última equação podemos resolvê-la usando fator integrante: ��

�� ,

�� ,

�� ,

Integrando dos dois lados:

��

�� ,

Logo, a solução geral de (IV) é:

�� 1� � �

�� Voltando em (II) chegamos à solução geral de (I):

��

Para determinar o lim�#$ �� aplicamos L´Hospital e obtemos: lim�#$�.��&'()

��&'() � �. Logo :

lim�#$ �� .

2

Esse limite é aceitável, pois considerando que � é o tamanho da população e dado que todos os indivíduos saudáveis “*��" dessa população estão suscetíveis à infecção, concluímos que, quando � # ∞, os indivíduos saudáveis tenderão à passar para a situação ��. Ou seja, toda a população � será igual a ��.

QUESTÃO 2

CONSIDERAÇÕES INICIAIS

• Os códigos computacionais completos dos métodos a seguir encontram-se em anexo.

• Teremos na questão 3 a validação de cada método, onde será confirma a ordem de consistência.

MÉTODO DE EULER O método de Euler é um método tipo RK11 (Runge-Kutta de 1ª. Ordem e 1 estágio) explícito, o método é dado por:

Dados �0 � �., �0 � �., / � /., 0 � 0, 1 � 30, implementamos em C utilizando um “for” ou “while”: /*cálculo de h*/ h=(b-a)/n; /*método de euler*/ for (i=0;i<=n;i++){ t[i+1]=t[i]+h; F[i]=Y(t[i]); //aqui temos a solução exata para a verificação do código k1= f(t[i],y[i]);//função a ser avaliada y[i+1]=y[i]+h*k1; }

MÉTODO DE EULER APRIMORADO O método de Euler Aprimorado é um método tipo RK22 (Runge-Kutta de 2ª. Ordem e 2 estágios) explícito, o método é dado por:

3

Implementando em C, temos: for (i=0;i<=n;i++){ t[i+1]=t[i]+h; F[i]=Y(t[i]); k1= f(t[i],y[i]); k2= f(t[i]+h,y[i]+h*k1); y[i+1]=y[i]+h*0.5*(k1+k2); }

MÉTODO DE EULER IMPLÍCITO (TRAPÉZIO) O método de Euler Implícito é dado por:

Nesse caso temos algumas modificações, pois o método é implícito: temos que utilizar um algoritmo para calcular a solução da função implícita, aqui foi utilizado o método de Newton e utilizamos o método de Euler para aproximar a solução ��3�, sendo o ponto inicial do método de Newton. Implementando em C, temos: /*método de Euler implicito*/ for (i=0;i<=n;i++){ t[i+1]=t[i]+h; F[i]=Y(t[i]); /*metodo de Euler para calcular o ponto inicial do metodo de Newton*/ k1= f_euler(t[i],z[i]); z[i+1]=z[i]+h*k1; y[i+1]=newton( z[i], y[i], h, 0.0001, MAX); }

MÉTODO DE RUNGE-KUTTA-FEHLBERG (RKF) O método RKF utiliza os métodos RK45 e RK56, a idéia desse método é obter uma estimativa do erro local de discretização e com isso efetuar o controle do erro, bem como uma variação do passo (h). Esse método tem ordem de consistência 4.

MÉTODO ADAMS-BASHFORTH DE 4 PASSOS O método de Adams-Bashforth de 4 Passos é um método de passo múltiplo explícito de quarta ordem:

4

Como o método é de quatro passos devemos ter quatro valores anteriores conhecidos, onde um pode ser dado pela condição inicial e os outros três podem ser calculados por outro método de quarta ordem. Neste exercício computacional foi utilizado o RK44 (clássico). Em C temos: /*método RK44 (classico)- calculo de y[1], y[2] e y[3]*/ for (i=0;i<=3;i++){ t[i+1]=t[i]+h; F[i]=Y(t[i]); k1= f(t[i],y[i]); k2= f(t[i]+0.5*h,y[i]+0.5*h*k1); k3= f(t[i]+0.5*h,y[i]+0.5*h*k2); k4= f(t[i]+h,y[i]+h*k3); y[i+1]=y[i]+h*(k1+2*k2+2*k3+k4)/6.0; } /* Adams-bashforth */ for (i=3;i<=n;i++){ t[i+1]=t[i]+h; F[i]=Y(t[i]); y[i+1]=y[i]+h*(55*f(t[i],y[i])-59*f(t[i-1],y[i-1])+37*f(t[i-2],y[i-2])-9*f(t[i-3],y[i-3]))/24; }

MÉTODO DE ADAMS-MOULTON DE 4 PASSOS O método de Adams-Moulton de 4 Passos é um método de passo múltiplo implícito de quinta ordem:

Para o método de Adams-Moulton de 4 passos foi utilizado um método RK56 para os cálculos dos �4 pré-

determinados. E como é um método implícito, foi necessário utilizar o método de Newton para se obter a raiz da equação obtida após a discretização. Essa raiz nos leva ao resultado de ��3�. Em C temos: /*método RK56 - calculo de y[1], y[2] e y[3]*/ for (i=0;i<=3;i++){ t[i+1]=t[i]+h; F[i]=Y(t[i]); k1 = f(t[i] , y[i]); k2 = f(t[i] + h/4 , y[i] + h*k1/4); k3 = f(t[i] + 3*h/8, y[i] + 3*h*k1/32 + 9*h*k2/32); k4 = f(t[i] + 12*h/13 , y[i] + 1932*h*k1/2197 - 7200*h*k2/2197 + 7296*h*k3/2197 ); k5 = f(t[i] + h , y[i] + 439*h*k1/216 - 8*h*k2 + 3680*h*k3/513 - 845*h*k4/4104); k6 = f(t[i] + h/2 , y[i] - 8*h*k1/27 + 2*h*k2 - 3544*h*k3/2565 + 1859*h*k4/4104 - 11*h*k5/40); y[i+1]=y[i]+h*(16*k1/135 + 6656*k3/12825 + 28561*k4/56430 - 9*k5/50 + 2*k6/55); } /* Adams-Moulton */ for (i=3;i<=n;i++){ t[i+1]=t[i]+h;

5

F[i]=Y(t[i]); /*metodo de Euler para calcular o ponto inicial do metodo de Newton*/ w1= f_euler(t[i],z[i]); z[i+1]=z[i]+h*w1; y[i+1]=newton( z[i], y[i], y[i-1], y[i-2] , y[i-3] , h, 0.00000000001, MAX); }

QUESTÃO 3


• Foi utilizada a solução exata de um Problema de Cauchy para fazer uma verificação no sentido de saber se o código foi implementado corretamente. Essa verificação ou também validação é baseada na análise da consistência do método. Esse procedimento será efetuado em todas as implementações exceto no método de RKF, tendo em vista que o mesmo possui controle automático de passo. O procedimento consiste na seguinte estratégia: Calculamos a relação:

Onde 5��; 7� é a solução aproximada pelo método e ��é a solução exata; o módulo |5��; 7� � ��| é o erro global no instante �.

Onde 9 é a ordem de consistência do método. Em C podemos ter o seguinte: val[j]=fabs(F[i-1]-y[i-1]); printf("Ordem de Consistencia p[%d]=%10.12f\n",j, log2(val[j-1]/val[j]));

• A validação foi feita utilizando a equação que foi dada no exercício computacional (equação 3):

��

Em C a função :��, �� é dada por: double f(double t, double y){ double Y; Y=k*(m-y)*y; return Y; }

E a solução exata conforme calculada na questão 1, em C fica da seguinte forma: double Y(double t){

6

double y; y=m*exp(k*m*t)/(exp(k*m*t)+m*C); return y; }

As constants são: double m=pow(10,5); double k=2*pow(10,-6); double C=99*pow(10,-5); double t0=0; double y0=pow(10,3);

A constante C é calculada utilizando a solução exata:

��

e impondo as condições iniciais: �., �.. Nesse caso obtemos � � 99.10�<.

• As tabelas de validação mostradas aqui, foram construídas por meio do Microsoft Excel com os dados obtidos dos resultados de saídas gerados pelos algoritmos implementados em C.

ITEM A

MÉTODO DE EULER Segue abaixo os resultados obtidos na comprovação numérica da ordem de consistência do método, ou seja, a validação:

h i=n ti sol. aprox sol. exata erro global consistência

0,6 50 30 77452,728 80295,71528 2842,987264

0,3 100 30 78928,111 80295,71528 1367,603849 1,1177481

0,15 200 30 79625,069 80295,71528 670,6459911 1,0615615

0,075 400 30 79963,637 80295,71528 332,0786715 1,0313549

0,0375 800 30 80130,481 80295,71528 165,2341043 1,0158094

0,01875 1600 30 80213,299 80295,71528 82,41649063 1,0079363

Com a tabela anterior observamos que a comprovação numérica da ordem de consistência tende a 1 conforme 7 diminui.

MÉTODO DE EULER APRIMORADO Na validação obtemos a seguinte tabela:

h i=n Ti sol. aprox sol. exata erro global consistência

0,6 50 30 80165,89826 80295,71528 129,8170159

0,3 100 30 80262,34633 80295,71528 33,36894391 1,994993

0,15 200 30 80287,25032 80295,71528 8,46495774 1,997712

0,075 400 30 80293,58313 80295,71528 2,132151677 1,998879

0,0375 800 30 80295,18021 80295,71528 0,535064302 1,999441

0,01875 1600 30 80295,58126 80295,71528 0,134021762 1,99972

7

Comprovamos que a ordem de consistência tende a 2.

MÉTODO DE EULER IMPLÍCITO Como explicado na questão 2, para aproximar a solução da equação (3), temos abaixo a implementação em C do método de Newton utilizando o método de Euler para aproximar a solução ��3� no método de Newton: /*Função a ser avaliada no EULER para o cálculo do ponto inicial no método de Newton*/ double f_euler(double tz, double z){ double Y; Y=k*(m-z)*z; return Y; } /*FUNCAO IMPLICITA */ double f(double x0, double yk, double h){ double x=0; x=-(h*k*pow(x0,2)*0.5) +(h*k*m*0.5-1)*x0-(h*k*pow(yk,2)*0.5)+(h*k*m*yk*0.5)+yk; return (x); } /*DERIVADA DA FUNCAO IMPLICITA*/ double Df(double x0, double h){ double dx; dx=-h*k*x0+h*k*m*0.5-1; return (dx); } /*METODO DE NEWTON*/ double newton(double x0, double yk, double h, double tol, int N){ int i; double x; for(i=0;1<=N;i++){ x=x0-f(x0,yk,h)/Df(x0,h); if(fabs(x-x0)<tol) return (x); else x0=x; } printf("o metodo falhou apos N="); return (N); }

Fazendo a validação do método obtemos:


0,6 50 30 80319,307 80295,71528 23,5917249

0,3 100 30 80301,61205 80295,71528 5,896770541 2,050244

0,15 200 30 80297,1894 80295,71528 1,47412021 2,024801

0,075 400 30 80296,0838 80295,71528 0,368525528 2,012316

0,0375 800 30 80295,80741 80295,71528 0,092131098 2,006136

0,01875 1600 30 80295,73831 80295,71528 0,023032755 2,003062

Verificamos que a ordem é 2.

MÉTODO DE RUNGE-KUTTA-FEHLBERG (RKF)

8

No método RKF não convém tentar efetuar a validação, pois como o método faz um controle do erro de discretização local, conseqüentemente, controlando também o erro global, temos uma variação do passo (h), bem como, variação no valor de (n). Podemos verificar esse fato na simulação abaixo: n=100 h inicial=0,30 i ti sol. aprox sol. exata erro global i ti sol. aprox sol. exata erro global

0 0 1000 1000 0 33 13,8 13762,95162 13762,94823 0,003392854

1 0,3 1061,18035 1061,180349 8,64161E-07 34 14,1 14490,68158 14490,67804 0,003546228

2 0,9 1194,860918 1194,86089 2,80979E-05 35 14,4 15250,08633 15250,08263 0,003702984

3 1,5 1345,152734 1345,152672 6,18202E-05 36 14,7 16041,82259 16041,81872 0,003862815

4 2,1 1514,058854 1514,058751 0,000103118 37 15 16866,4829 16866,47887 0,004025374

5 2,7 1703,807637 1703,807484 0,000153214 38 15,3 17724,58671 17724,58252 0,004190269

6 3,3 1916,873761 1916,873547 0,000213482 39 15,6 18616,57125 18616,5669 0,004357061

7 3,9 2156,000012 2155,999727 0,000285444 40 15,9 19542,78226 19542,77773 0,004525268

8 4,5 2424,219504 2424,219133 0,000370783 41 16,2 20503,46463 20503,45994 0,004694362

9 5,1 2724,877852 2724,877381 0,00047134 42 16,5 21498,75325 21498,74838 0,004863772

10 5,7 3061,654664 3061,654075 0,000589107 43 16,8 22528,6639 22528,65886 0,005032882

11 6,3 3438,583505 3438,582779 0,000726218 44 17,1 23593,08465 23593,07944 0,00520104

12 6,9 3860,069275 3860,06839 0,000884925 45 17,7 25824,32177 25824,31617 0,005606576

13 7,5 4330,901646 4330,900578 0,001067568 46 18,3 28188,72232 28188,71633 0,005992486

14 8,1 4856,262923 4856,261646 0,00127653 47 18,9 30680,06246 30680,05611 0,00635083

15 8,4 5141,118955 5141,117605 0,001350473 48 19,5 33289,51733 33289,51065 0,006673935

16 8,7 5441,728254 5441,726826 0,001428126 49 20,1 36005,6357 36005,62874 0,006954713

17 9 5758,847551 5758,846041 0,001509609 50 20,7 38814,42395 38814,41676 0,007186981

18 9,3 6093,256272 6093,254677 0,001595033 51 21,3 41699,54824 41699,54087 0,007365769

19 9,6 6445,755501 6445,753817 0,001684507 52 22,5 47623,8027 47623,79509 0,00760509

20 9,9 6817,16668 6817,164901 0,00177813 53 23,1 50621,97927 50621,97166 0,007609775

21 10,2 7208,33001 7208,328134 0,001875992 54 23,7 53615,68901 53615,68145 0,007557599

22 10,5 7620,102537 7620,100559 0,001978172 55 24,3 56583,56346 56583,55601 0,007452555

23 10,8 8053,355887 8053,353802 0,002084738 56 24,9 59504,96674 59504,95944 0,007300337

24 11,1 8508,973623 8508,971428 0,002195741 57 25,5 62360,55684 62360,54973 0,007107998

25 11,4 8987,84821 8987,845899 0,002311217 58 26,1 65132,7716 65132,76472 0,006883501

26 11,7 9490,877547 9490,875116 0,002431181 59 26,7 67806,21539 67806,20876 0,00663519

27 12 10018,96105 10018,9585 0,002555628 60 27,3 70367,93118 70367,92481 0,006371255

28 12,3 10572,99529 10572,99261 0,002684529 61 27,9 72807,55318 72807,54708 0,006099243

29 12,6 11153,86909 11153,86627 0,002817829 62 28,5 75117,34397 75117,33815 0,005825675

30 12,9 11762,4582 11762,45524 0,002955444 63 29,1 77292,12803 77292,12247 0,005555794

31 13,2 12399,61939 12399,61629 0,003097258 64 29,7 79329,13877 79329,13348 0,005293464

32 13,5 13066,18412 13066,18087 0,003243123 65 30 80295,72039 80295,71528 0,005110912

Na tabela anterior podemos ver que houve uma variação no valor de (n), foi dado como entrada (n=100) e tivemos na saída (n=65). Essa variação de (n) ocorre devido a variação de (h), colocando a variação de (h) em evidência em cada instante de mudança, temos:

i h

0 0,3

1 0,6

14 0,3

44 0,6

51 1,2

52 0,6

57 1,2

65 0,3

9

MÉTODO ADAMS-BASHFORTH DE 4 PASSOS Abaixo temos a validação do método


0,6 50 30 80295,41941 80295,71528 0,295870682

0,3 100 30 80295,69984 80295,71528 0,015433206 4,384676

0,15 200 30 80295,71439 80295,71528 0,000882245 4,251047

0,075 400 30 80295,71522 80295,71528 5,28443E-05 4,141201

0,0375 800 30 80295,71527 80295,71528 3,23728E-06 4,074058

0,01875 1600 30 80295,71528 80295,71528 2,02679E-07 4,021792

Com isso verificamos que de fato a ordem de consistência do método é 4.

MÉTODO DE ADAMS-MOULTON DE 4 PASSOS As funções que são utilizadas no método de Newton para o método de Adams-Moulton de 4 passos são similares as do método de Euler Implícito, podendo ser encontradas no anexo desse trabalho. No decorrer da implementação desse método, vale ressaltar que, o método se comportou de forma estranha no começo das simulações (valores completamente anormais na consistência e não havia diminuição do erro global), fazendo alguns testes foi percebido que era necessário diminuir ainda mais a tolerância do erro do método de Newton chegando ao valor (0,00000000001):

y[i+1]=newton( z[i], y[i], y[i-1], y[i-2] , y[i-3] , h, 0.00000000001, MAX);. Dessa forma, foi conseguido resultados um pouco melhores nas simulações, como se pode observar na validação:


3 10 30 80285,10144 80295,71528 10,61383366

1,5 20 30 80294,74005 80295,71528 0,975222323 4,118352

0,75 40 30 80295,67949 80295,71528 0,035791832 4,272095

0,375 80 30 80295,71413 80295,71528 0,001146339 4,903719

0,1875 160 30 80295,71525 80295,71528 3,05982E-05 5,224729

0,09375 320 30 80295,71528 80295,71528 4,41719E-06 2,808191

Na tabela acima vemos que os resultados obtidos na coluna “consistência” se comportam bem até / � 160, depois disso passamos a ter valores estranhos para o teste da consistência. Entretanto, observamos que o erro global continua a diminuir conforme / aumenta, predizendo a boa estabilidade do método.

OTIMIZANDO O PASSO TEMPORAL (h) Para fazer uma otimização do passo temporal neste trabalho, primeiramente foi estabelecido um critério. Esse critério foi relacionado com a precisão da solução aproximada por cada método, tomando como base o erro global que foi calculado. A fim de se ter uma idéia da magnitude da precisão requerida, foi analisado o problema proposto pelo Exercício Computacional. Verificando o significado da solução ��, a mesma indica a quantidade de indivíduos infectados ao longo do tempo, logo parece razoável tomar como precisão a primeira casa decimal do erro global, tendo em vista que a contagem dos indivíduos se dá no conjunto dos números inteiros. Analisando rapidamente as tabelas de validação obtidas anteriormente, podemos ter uma idéia inicial do passo (h) que deve ser tomado em cada método:

10

Método h n Erro global

EULER 0,01875 1600 82,41649063

EULER APRIMORADO 0,0375 800 0,535064302

EULER IMPLÍCITO 0,075 400 0,368525528

ADAMS-BASHOTH DE 4 PASSOS 0,6 50 0,295870682

ADMS-MOULTON DE 4 PASSOS 1,5 20 0,975222323

Essa tabela nos dá uma idéia do tamanho do passo que devemos tomar levando em conta a precisão requerida. No entanto, observando o método de Euler vemos que o erro global não atinge a precisão requerida, isso ocorre pelo fato do método de Euler ser de consistência 1, onde seria necessário um (n) muito grande (h muito pequeno) para ser atingida a casa decimal do erro global. Ressalta-se que na tabela acima temos apenas uma base de um (h) inicial a ser adotado, para se ter um (h) mais apurado seria interessante implementar um algoritmo para o controle do passo de integração. Esse procedimento é feito no método de Runge-Kutta-Fehlberg (RKF).

ITEM B


• As tabelas mostradas a seguir foram geradas pelos códigos computacionais no anexo, para a implementação dos códigos foi utilizado o programa Dev-C++ versão 4.9.9.2.

• Os gráficos foram plotados utilizando os dados das tabelas através do programa "/>?@�".

• Nos gráficos temos a seguinte legenda:

MÉTODO DE EULER SIMULAÇÃO 1 Número de iterações

50

h=0,6

i ti sol. aprox sol. exata erro global i ti sol. aprox sol. exata erro global

0 0 1000 1000 0 26 15,6 16377,49046 18616,5669 2239,076439

1 0,6 1118,8 1126,06116 7,2611591 27 16,2 18020,92268 20503,45994 2482,537258

2 1,2 1251,55394 1267,81023 16,256282 28 16,8 19793,72902 22528,65886 2734,929849

3 1,8 1399,86075 1427,14517 27,284422 29 17,4 21698,82645 24691,76229 2992,935837

4 2,4 1565,49251 1606,17914 40,686633 30 18 23737,67874 26990,19983 3252,521092

5 3 1750,41069 1807,26097 56,850277 31 18,6 25910,02732 29419,00856 3508,981242

6 3,6 1956,78325 2032,99662 76,21337 32 19,2 28213,63518 31970,6736 3757,038425

7 4,2 2187,00244 2286,27126 99,268818 33 19,8 30644,06035 34635,05355 3990,993205

8 4,8 2443,70316 2570,2715 126,56834 34 20,4 33194,47747 37399,40822 4204,930754

9 5,4 2729,78151 2888,50734 158,72583 35 21 35855,56676 40248,5411 4392,974336

10 6 3048,41325 3244,833 196,41976 36 21,6 38615,48877 43165,06136 4549,572594

11 6,6 3403,07145 3643,46576 240,39431 37 22,2 41459,96026 46129,76116 4669,800908

12 7,2 3797,54295 4089,00152 291,45857 38 22,8 44372,44152 49122,09398 4749,652458

13 7,8 4235,9425 4586,4257 350,4832 39 23,4 47334,43822 52120,73064 4786,292417

11

14 8,4 4722,72375 5141,1176 418,39385 40 24 50325,91196 55104,16192 4778,249959

15 9 5262,68566 5758,84604 496,16038 41 24,6 53325,7845 58051,31195 4725,52745

16 9,6 5860,9729 6445,75382 584,78091 42 25,2 56312,51149 60942,12564 4629,614157

17 10 6523,06845 7208,32813 685,25969 43 25,8 59264,69413 63758,09591 4493,401788

18 11 7254,77615 8053,3538 798,57765 44 26,4 62161,69266 66482,70245 4321,009791

19 11 8062,19116 8987,8459 925,65474 45 27 64984,20453 69101,74235 4117,537817

20 12 8951,65539 10018,9585 1067,3031 46 27,6 67714,77287 71603,54255 3888,769675

21 13 9929,69547 11153,8663 1224,1708 47 28,2 70338,19706 73979,05361 3640,856554

22 13 11002,9403 12399,6163 1396,676 48 28,8 72841,82635 76221,83323 3380,006881

23 14 12178,0155 13762,9482 1584,9327 49 29,4 75215,72751 78327,93408 3112,206568

24 14 13461,4125 15250,0826 1788,6701 50 30 77452,72801 80295,71528 2842,987264

25 15 14859,3304 16866,4789 2007,1484

TABELA 1 Nessa tabela, embora o erro global passe a diminuir após A � 44, ilustrando uma estabilidade, observamos um erro global muito grande. Podemos observar claramente no gráfico 1 a ilustrações do erro global que é cometido entre a solução exata e solução aproximada.

SIMULAÇÃO 2

Número de iteraçoes

1600

h=0,018750

i ti sol. aprox sol. exata erro global i ti sol, aprox sol, exata erro global

0 0 1000 1000 0 840 15,75 18996,12924 19075,37623 79,246991

40 0,75 1159,6423 1159,957035 0,31469485 880 16,5 21409,27881 21498,74838 89,4695752

80 1,5 1344,4255 1345,152672 0,72713631 920 17,25 24038,26838 24138,15909 99,8907128

120 2,25 1558,1905 1559,449643 1,25910783 960 18 26879,97579 26990,19983 110,224045

160 3 1805,3247 1807,260968 1,9362746 1000 18,75 29925,62251 30045,75813 120,135623

200 3,75 2090,8257 2093,614344 2,78865136 1040 19,5 33160,25109 33289,51065 129,259563

240 4,5 2420,3681 2424,219133 3,85104837 1080 20,25 36562,51757 36699,73805 137,22048

280 5,25 2800,3706 2805,534083 5,16345571 1120 21 40104,87994 40248,5411 143,661157

320 6 3238,0617 3244,833003 6,77131048 1160 21,75 43754,23067 43902,50325 148,272584

360 6,75 3741,5389 3750,264515 8,72557154 1200 22,5 47472,97256 47623,79509 150,822537

400 7,5 4319,8181 4330,900578 11,0825039 1240 23,25 51220,48099 51371,6595 151,178509

440 8,25 4982,8638 4996,766854 13,903049 1280 24 54954,84058 55104,16192 149,321336

480 9 5741,5944 5758,846041 17,2516319 1320 24,75 58634,70278 58780,05002 145,347236

12

520 9,75 6607,849 6629,043272 21,1942354 1360 25,5 62221,09183 62360,54973 139,4579

560 10,5 7594,305 7620,100559 25,7955596 1400 26,25 65678,99459 65810,93491 131,940323

600 11,25 8714,3304 8745,445523 31,1150924 1440 27 68978,60273 69101,74235 123,139624

640 12 9981,7565 10018,9585 37,2019598 1480 27,75 72096,12658 72209,55557 113,428989

680 12,75 11410,554 11454,64228 44,0885068 1520 28,5 75014,15737 75117,33815 103,180775

720 13,5 13014,398 13066,18087 51,7827065 1560 29,25 77721,60748 77814,34951 92,742039

760 14,25 14806,119 14866,37854 60,2597101 1600 30 80213,29879 80295,71528 82,4164906

800 15 16797,026 16866,47887 69,4531313

TABELA 2 Na tabela 2 observamos a diminuição do erro global após a diminuição de (h), que também pode ser observado no gráfico:

No gráfico 2, devido à escala temos a impressão que não existe diferença entre solução aproximada e solução exata. Essa diferença pode ser melhor observada na tabela 2.

13

MÉTODO DE EULER APRIMORADO

Comparando o gráfico 1 com o Gráfico 3 observamos uma semelhança, no entanto no gráfico 1 temos o método de Euler com / � 50 e no gráfico 2 o método de Euler Aprimorado com / � 10, com isso podemos constatar a vantagem do método e Euler Aprimorado em relação ao Euler. SIMULAÇÃO 2

Número de iterações 1600

h=0,018750


0 0 1000 1000 0 840 15,75 19075,28 19075,38 0,0940974

40 0,8 1159,957 1159,957 0,00038909 880 16,5 21498,64 21498,75 0,1059267

80 1,5 1345,152 1345,153 0,00089838 920 17,25 24138,04 24138,16 0,1179656

120 2,3 1559,448 1559,45 0,00155437 960 18 26990,07 26990,2 0,1299151

160 3 1807,259 1807,261 0,00238817 1000 18,75 30045,62 30045,76 0,1414385

200 3,8 2093,611 2093,614 0,00343598 1040 19,5 33289,36 33289,51 0,1521795

240 4,5 2424,214 2424,219 0,00473961 1080 20,25 36699,58 36699,74 0,1617857

280 5,3 2805,528 2805,534 0,0063468 1120 21 40248,37 40248,54 0,1699339

320 6 3244,825 3244,833 0,00831145 1160 21,75 43902,33 43902,5 0,176355

360 6,8 3750,254 3750,265 0,01069354 1200 22,5 47623,61 47623,8 0,1808546

400 7,5 4330,887 4330,901 0,01355868 1240 23,25 51371,48 51371,66 0,1833269

440 8,3 4996,75 4996,767 0,01697719 1280 24 55103,98 55104,16 0,1837587

480 9 5758,825 5758,846 0,0210224 1320 24,75 58779,87 58780,05 0,1822253

520 9,8 6629,018 6629,043 0,02576814 1360 25,5 62360,37 62360,55 0,1788778

560 11 7620,069 7620,101 0,03128523 1400 26,25 65810,76 65810,93 0,1739251

600 11 8745,408 8745,446 0,03763669 1440 27 69101,57 69101,74 0,1676142

640 12 10018,91 10018,96 0,04487175 1480 27,75 72209,4 72209,56 0,1602101

14

680 13 11454,59 11454,64 0,05301871 1520 28,5 75117,19 75117,34 0,1519783

720 14 13066,12 13066,18 0,06207676 1560 29,25 77814,21 77814,35 0,143172

760 14 14866,31 14866,38 0,07200748 1600 30 80295,58 80295,72 0,1340218

800 15 16866,4 16866,48 0,08272661

TABELA 4 Comparando mais uma vez com o Euler temos uma maior precisão no erro global para o método de Euler Aprimorado. Graficamente temos:

EULER IMPLÍCITO (TRAPÉZIO) Simulação 1


h=3,00

i ti sol, aprox sol, exata erro global

0 0 1000 1000 0

1 3 1838,37308 1807,261 31,1121079

2 6 3351,49791 3244,833 106,664908

3 9 6020,71791 5758,846 261,871868

4 12 10549,055 10018,958 530,096457

5 15 17762,0717 16866,479 895,592826

6 18 28221,282 26990,2 1231,08212

7 21 41585,9559 40248,541 1337,41476

8 24 56256,1492 55104,162 1151,98731

9 27 69945,2875 69101,742 843,545186

10 30 80889,3822 80295,715 593,666969

TABELA 5 Comparando a tabela 5 com a tabela 3 notamos a diminuição do erro global no método de Euler implicito, caracterizando a vantagem desse método.

15

Simulação 2 Número de iterações 50 H=0,60

i ti sol, aprox sol, exata erro global i ti sol, aprox sol, exata erro global

0 0 1000 1000 0 26 15,6 18654,128 18616,57 37,560903

1 0,6 1126,21169 1126,0612 0,15053067 27 16,2 20543,996 20503,46 40,535932

2 1,2 1268,1473 1267,8102 0,33707297 28 16,8 22572,026 22528,66 43,367039

3 1,8 1427,71085 1427,1452 0,56567979 29 17,4 24737,741 24691,76 45,978311

4 2,4 1607,02232 1606,1791 0,84317379 30 18 27038,493 26990,2 48,293234

5 3 1808,43816 1807,261 1,17719656 31 18,6 29469,248 29419,01 50,238947

6 3,6 2034,57287 2032,9966 1,57624625 32 19,2 32022,424 31970,67 51,750894

7 4,2 2288,32095 2286,2713 2,0496979 33 19,8 34687,831 34635,05 52,777501

8 4,8 2572,8793 2570,2715 2,60779927 34 20,4 37452,693 37399,41 53,284367

9 5,4 2891,76897 2888,5073 3,2616329 35 21 40301,799 40248,54 53,257469

10 6 3248,85604 3244,833 4,02303384 36 21,6 43217,766 43165,06 52,704919

11 6,6 3648,37021 3643,4658 4,90445053 37 22,2 46181,418 46129,76 51,656962

12 7,2 4094,92025 4089,0015 5,91873446 38 22,8 49172,258 49122,09 50,164112

13 7,8 4593,50454 4586,4257 7,0788434 39 23,4 52169,024 52120,73 48,293557

14 8,4 5149,51505 5141,1176 8,3974423 40 24 55150,286 55104,16 46,124211

15 9 5768,73243 5758,846 9,8863866 41 24,6 58095,053 58051,31 43,740984

16 9,6 6457,30989 6445,7538 11,5560756 42 25,2 60983,355 60942,13 41,228913

17 10,2 7221,7428 7208,3281 13,4146681 43 25,8 63796,764 63758,1 38,667784

18 10,8 8068,82096 8053,3538 15,467161 44 26,4 66518,83 66482,7 36,127785

19 11,4 9005,56024 8987,8459 17,7143441 45 27 69135,409 69101,74 33,666506

20 12 10039,1102 10018,958 20,1516598 46 27,6 71634,87 71603,54 31,327428

21 12,6 11176,6343 11153,866 22,7680189 47 28,2 74008,193 73979,05 29,139794

22 13,2 12425,1609 12399,616 25,5446475 48 28,8 76248,953 76221,83 27,119651

23 13,8 13791,4023 13762,948 28,4540678 49 29,4 78353,206 78327,93 25,271717

24 14,4 15281,542 15250,083 31,4593413 50 30 80319,307 80295,72 23,591725

25 15 16900,9926 16866,479 34,5137258

TABELA 6

16

MÉTODO DE ADAMS-BASHFORTH DE 4 PASSOS


h=0,60 i ti sol. aprox sol. exata erro global i ti sol. aprox sol. exata erro global

0 0 1000 1000 0 26 16 18616,876 18616,567 0,30935601

1 0,6 1126,061 1126,061 0,0001984 27 16 20503,899 20503,46 0,4388509

2 1,2 1267,81 1267,81 0,00044465 28 17 22529,236 22528,659 0,57669505

3 1,8 1427,144 1427,145 0,00074694 29 17 24692,479 24691,762 0,71626508

4 2,4 1606,171 1606,179 0,00781103 30 18 26991,049 26990,2 0,84923582

5 3 1807,244 1807,261 0,01706095 31 19 29419,975 29419,009 0,96601872

6 3,6 2032,97 2032,997 0,02705828 32 19 31971,73 31970,674 1,05649433

7 4,2 2286,233 2286,271 0,03846274 33 20 34636,165 34635,054 1,11100946

8 4,8 2570,22 2570,271 0,05134955 34 20 37400,53 37399,408 1,12154836

9 5,4 2888,442 2888,507 0,06550691 35 21 40249,624 40248,541 1,08292619

10 6 3244,752 3244,833 0,08075246 36 22 43166,055 43165,061 0,99380859

11 6,6 3643,369 3643,466 0,09681268 37 22 46130,619 46129,761 0,85735201

12 7,2 4088,888 4089,002 0,11325221 38 23 49122,775 49122,094 0,68129566

13 7,8 4586,296 4586,426 0,12946347 39 23 52121,208 52120,731 0,47741915

14 8,4 5140,973 5141,118 0,14464137 40 24 55104,422 55104,162 0,26039603

15 9 5758,688 5758,846 0,15775319 41 25 58051,358 58051,312 0,04619531

16 9,6 6445,586 6445,754 0,16751819 42 25 60941,976 60942,126 0,14972019

17 10 7208,156 7208,328 0,17240001 43 26 63757,782 63758,096 0,3141059

18 11 8053,183 8053,354 0,17061736 44 26 66482,265 66482,702 0,43728576

19 11 8987,686 8987,846 0,16018158 45 27 69101,228 69101,742 0,51394814

20 12 10018,82 10018,96 0,13897035 46 28 71602,999 71603,543 0,54331014

21 13 11153,76 11153,87 0,10484632 47 28 73978,525 73979,054 0,52869362

22 13 12399,56 12399,62 0,05582798 48 29 76221,357 76221,833 0,47665305

23 14 13762,96 13762,95 0,00968399 49 29 78327,538 78327,934 0,39584928

24 14 15250,18 15250,08 0,09262837 50 30 80295,419 80295,715 0,29587068

25 15 16866,67 16866,48 0,1929668

TABELA 7

17

Observando as tabelas 6 E 7 Fica clara a melhor precisão do método de Adams-Bashforth em relação ao Euler Implícto.

MÉTODO RUNGE-KUTTA-FEHLBERG SIMULAÇÃO 1

Número de passos n=100

TABELA 8 A tabela 8 nos traz alguns pontos de uma simulação do método RKF com / � 100. Observamos que devido ao controle do passo, na saída temos / � 65, Como o método RKF é de quarta ordem temos uma ótima precisão mostrada no baixo valor do erro global. Fazendo a comparação com o Euler implícito o RKF é extremamente vantajoso e melhor também que o método de Adams-Bashforth de 4 passos. Plotando os pontos acima temos:

MÉTODO DE ADAMS-MOULTON DE 4 PASSOS

Número de iterações

50

18

h=0,60


0 0 1000 1000 0 26 15,6 18616,556 18616,567 0,0108427

1 0,6 1126,0612 1126,0612 2,55732E-06 27 16,2 20503,448 20503,46 0,0117715

2 1,2 1267,8102 1267,8102 5,72372E-06 28 16,8 22528,646 22528,659 0,0124187

3 1,8 1427,1452 1427,1452 9,60082E-06 29 17,4 24691,75 24691,762 0,0126922

4 2,4 1606,1792 1606,1791 9,89027E-06 30 18 26990,187 26990,2 0,0125162

5 3 1807,261 1807,261 2,83207E-05 31 18,6 29418,997 29419,009 0,011846

6 3,6 2032,9967 2032,9966 4,24835E-05 32 19,2 31970,663 31970,674 0,0106811

7 4,2 2286,2713 2286,2713 4,97181E-05 33 19,8 34635,044 34635,054 0,0090761

8 4,8 2570,2715 2570,2715 4,60323E-05 34 20,4 37399,401 37399,408 0,0071443

9 5,4 2888,5074 2888,5073 2,6358E-05 35 21 40248,536 40248,541 0,005052

10 6 3244,833 3244,833 1,56423E-05 36 21,6 43165,058 43165,061 0,0030035

11 6,6 3643,4657 3643,4658 8,77388E-05 37 22,2 46129,76 46129,761 0,0012161

12 7,2 4089,0013 4089,0015 0,00019924 38 22,8 49122,094 49122,094 0,0001095

13 7,8 4586,4253 4586,4257 0,00036102 39 23,4 52120,731 52120,731 0,0008202

14 8,4 5141,117 5141,1176 0,000585423 40 24 55104,163 55104,162 0,0008345

15 9 5758,8452 5758,846 0,00088598 41 24,6 58051,312 58051,312 0,0001562

16 9,6 6445,7525 6445,7538 0,001276881 42 25,2 60942,125 60942,126 0,0011268

17 10,2 7208,3264 7208,3281 0,001772112 43 25,8 63758,093 63758,096 0,0028566

18 10,8 8053,3514 8053,3538 0,002384208 44 26,4 66482,698 66482,702 0,0048312

19 11,4 8987,8428 8987,8459 0,003122551 45 27 69101,736 69101,742 0,0068347

20 12 10018,955 10018,958 0,003991207 46 27,6 71603,534 71603,543 0,0086671

21 12,6 11153,861 11153,866 0,004986363 47 28,2 73979,043 73979,054 0,010168

22 13,2 12399,61 12399,616 0,006093516 48 28,8 76221,822 76221,833 0,011231

23 13,8 13762,941 13762,948 0,007284717 49 29,4 78327,922 78327,934 0,0118081

24 14,4 15250,074 15250,083 0,008516316 50 30 80295,703 80295,715 0,0119049

25 15 16866,469 16866,479 0,009727837

TABELA 9

19

Analisando a tabela 9 o método de Adams-Moulton se apresenta com melhor precisão que o Adams-Bashforth mas em relação ao RKF essa comparação torna-se difícil, pois embora o Adams-Moulton seja de quinta ordem o método RKF (de quarta ordem) possui o controle do passo, fato que também é bastante vantajoso. Em relação à estabilidade pode-se dizer que o método de Adams-Molton, bem como, os outros métodos, se comportaram bem, ou seja, o erro global sempre se manteve controlado, embora às vezes com valor elevado sendo que na maioria dos casos em dado instante o erro passava a diminuir.

ANÁLISE DO CUSTO COMPUTACIONAL DE CADA MÉTODO E CONCLUSÕES FINAIS Tendo em vista o problema proposto pelo Exercício Computacional e os critérios estabelecidos para a precisão de cada método (primeira casa decimal do erro global), utilizando o (h) de forma otimizada conforme a análise da questão 3 item (b), temos o seguinte resultado em relação ao tempo de execução de cada método (custo computacional medido pelo tempo de execução do programa):

Colocação Método

Ordem de consistencia h n

Erro global em t=n tempo

1º ADAMS-BASHFORTH DE 4 PASSOS 4 0,6 50 0,295870682 0.015000 s

2º RKF 4 VARIÁVEL 100 0.005110911792 0.031000 s

3º EULER IMPLÍCITO 2 0,075 400 0,368525528 0.297000 s

4º EULER APRIMORADO 2 0,0375 800 0,535064302 0.625000 s

5º ADAMS-MOULTON DE 4 PASSOS 5 1,5 20 0,975222323 1.453000 s

6º EULER 1 0,01875 1600 82,41649063 não avaliado

TABELA 10 O tempo em segundos foi obtido utilizando os seguintes comandos em C: clock_t to, tf; double tempo_gasto; to = clock(); . . . tf = clock(); tempo_gasto = ( (double) (tf - to) ) / CLOCKS_PER_SEC; printf("Tempo gasto: %lf s\n", tempo_gasto);

Na tabela 10 observamos que no método de Euler não foi avaliado o tempo de execução, pois nem mesmo fazendo / � 50000 (capacidade de iterações do programa) foi conseguido a precisão requerida. Em relação à colocação dos outros métodos vale ressaltar que para o problema proposto onde se tomou como critério de precisão a primeira casa decimal do erro global, acaba sendo mais vantajoso a utilização do método de Adams-Bashforth ou RKF os quais nos deram a precisão desejada no menor tempo de execução. Os método de Euler Implícito e Euler Aprimorado também tiveram bons desempenhos para o problema em questão. Por outro lado se tivéssemos um problema que requeresse uma maior precisão possivelmente o método de Adams-Moulton seria adequado às necessidades, pois ele ganha na ordem de consistência e conseqüentemente a precisão aumenta.

BIBLIOGRAFIA

1- A. M. Roma, R.L. Nós: Notas do Curso tratamento Numérico de Equações Diferenciais. Última revisão em 2011 – USP.

2- W. E. Boyce, R.C. Diprima: Equações Diferenciais Elementares e Problemas de Valores de Contorno”. 8ª. Ed. – LTC.

3- R. L. Burden, J. D. Faires: Análise Numérica. 7ª. Ed. – Thomson

20

ANEXOS – CÓDIGOS COMPUTACIONAIS /* Com exceção do método RKF (executa para n=100), os códigos abaixo requerem a entrada do número de simulações que serão realizadas (sim), e o valor de (n). Será fornecida uma saída em forma de tabela das simulações contendo as (n) iterações para t=b=30, caso necessite resultados de validações ou tempo de execução são necessárias algumas modificações.

Programa utilizado: Dev-C++ versão 4.9.9.2. */

/*METODO DE EULER*/ #include <stdio.h> #include<stdlib.h> #include <math.h> #define MAX 50000 double m=pow(10,5); double k=2*pow(10,-6); double C=99*pow(10,-5); double t0=0; double y0=pow(10,3); /*Função a ser avaliada*/ double f(double t, double y){ double Y; Y=k*(m-y)*y; return Y; } /*solucao exata*/ double Y(double t){ double y; y=m*exp(k*m*t)/(exp(k*m*t)+m*C); return y; } int main() { int i,j; //contadores int sim; //número de simulações int n; //quantidade de iterações double t[MAX]/*tempo*/, y[MAX]/*sol. aprox.*/, h /*tamanho do passo*/, k1, F[MAX]/*sol. exata*/, val[MAX], a=0, b=30; t[0]=t0;/*condição inicial*/ y[0]=y0;/*condição inicial*/ printf("Numero de simulacoes:\n"); scanf("%d", &sim); printf("\n"); for(j=0;j<sim;j++){ printf("Numero de iteracoes n:\n"); scanf("%d", &n); //printf("\n"); //printf("Valor de b - qtde de dias (30 dias)\n"); //scanf("%lf",&b); /*cálculo de h*/ h=(b-a)/n; printf("h=%f", h); printf("\n"); printf("\n"); /*método de euler*/ for (i=0;i<=n;i++){ t[i+1]=t[i]+h; F[i]=Y(t[i]); k1= f(t[i],y[i]);

21

y[i+1]=y[i]+h*k1; } // i=i-1; //usado na validacao /*fim do método*/ printf(" i ti sol. aprox sol. exata erro global\n"); printf("\n"); i=0; while(i<=n){ //while para a impressao da tabela com todos os i /*impressão da tabela*/ printf("%5d %10.3f %18.12f %18.12f %18.12f\n",i,t[i], y[i], F[i], fabs(F[i]-y[i])); i=i+1; } printf("\n"); printf("\n"); val[j]=fabs(F[i-1]-y[i-1]); printf("\n"); /*validação ->|sol. aprox(h)-sol. exata|/|sol. aprox(h/2)-sol. exata|*/ // printf("Ordem de Consistencia p[%d]=%10.12f\n",j, log2(val[j-1]/val[j-0])); // printf("\n"); } printf("\n"); system("pause"); return 0; }

/*EULER APRIMORADO*/ #include <stdio.h> #include<stdlib.h> #include <math.h> #include <time.h> #define MAX 50000 double m=pow(10,5); double k=2*pow(10,-6); double C=99*pow(10,-5); double t0=0; double y0=pow(10,3); /*Função*/ double f(double t, double y){ double Y; Y=k*(m-y)*y; return Y; } /*solucao exata*/ double Y(double t){ double y; y=m*exp(k*m*t)/(exp(k*m*t)+m*C); return y; } int main() { int i, n, j, sim; /*usado na verificação do tempo de execução clock_t to, tf; double tempo_gasto; to = clock();*/

22

double t[MAX], y[MAX]/*sol. aprox.*/, h, F[MAX]/*sol. exata*/, val[MAX], a, b=30, k1,k2; a=0; t[0]=t0;/*condição inicial*/ y[0]=y0;/*condição inicial*/ printf("Numero de simulacoes:\n"); scanf("%d", &sim); printf("\n"); for(j=0;j<sim;j++){ printf("Numero de iteracoes n:\n"); scanf("%d", &n); printf("\n"); //printf("Valor de b - qtde de dias (30 dias)\n"); //scanf("%lf",&b); h=(b-a)/n; printf("h=%f", h); printf("\n"); printf("\n"); /*método de Euler Aprimorado*/ for (i=0;i<=n;i++){ t[i+1]=t[i]+h; F[i]=Y(t[i]); k1= f(t[i],y[i]); k2= f(t[i]+h,y[i]+h*k1); y[i+1]=y[i]+h*0.5*(k1+k2); } // i=i-1;//utilizado para validacao do metodo /*fim do método*/ printf(" i ti sol. aprox sol. exata erro global\n"); printf("\n"); i=0; while(i<=n){ //for para a impressao da tabela com todos os valores de i /*impressão da tabela*/ printf("%5d %10.3f %18.12f %18.12f %18.12f\n",i,t[i], y[i], F[i], fabs(F[i]-y[i])); i=i+1; } printf("\n"); val[j]=fabs(F[i-1]-y[i-1]); printf("\n"); // /*validação ->|sol. aprox(h)-sol. exata|/|sol. aprox(h/2)-sol. exata|*/ // printf("Ordem de Consistencia p[%d]=%10.12f\n",j, log2(val[j-1]/val[j-0])); // printf("\n"); } /*usado na verificação do tempo tf = clock(); tempo_gasto = ( (double) (tf - to) ) / CLOCKS_PER_SEC; printf("Tempo gasto: %lf s\n", tempo_gasto);*/ printf("\n"); system("pause"); return 0; }

/*MÉTODO DE EULER IMPLÍCITO (TRAPÉZIO) #include <stdio.h> #include<stdlib.h> #include <math.h> #include <time.h> #define MAX 50000

23

double m=pow(10,5); double k=2*pow(10,-6); double C=99*pow(10,-5); double t0=0; double y0=pow(10,3); /*solucao exata*/ double Y(double t){ double y; y=m*exp(k*m*t)/(exp(k*m*t)+m*C); return y; } /*Função a ser avaliada no EULER para o caulo do ponto inicial no metodo de Newton*/ double f_euler(double tz, double z){ double Y; Y=k*(m-z)*z; return Y; } /*FUNCAO IMPLICITA*/ double f(double x0, double yk, double h){ double x=0; x=-(h*k*pow(x0,2)*0.5) +(h*k*m*0.5-1)*x0-(h*k*pow(yk,2)*0.5)+(h*k*m*yk*0.5)+yk; return (x); } /*DERIVADA DA FUNCAO IMPLICITA*/ double Df(double x0, double h){ double dx; dx=-h*k*x0+h*k*m*0.5-1; return (dx); } /*METODO DE NEWTON*/ double newton(double x0, double yk, double h, double tol, int N){ int i; double x; for(i=0;1<=N;i++){ x=x0-f(x0,yk,h)/Df(x0,h); if(fabs(x-x0)<tol) return (x); else x0=x; } printf("o metodo falhou apos N="); return (N); } int main() { int i=0, n, j, sim; /*usado para a verificação do tempo de execução clock_t to, tf; double tempo_gasto; to = clock();*/ double t[MAX], y[MAX]/*sol. aprox.*/, z[MAX], h, F[MAX]/*sol. exata*/, val[MAX], k1 , a, b=30; a=0; z[0]=y0;/*condicao inicial Euler*/ t[0]=t0;/*condição inicial*/ y[0]=y0;/*condição inicial*/

24

printf("Numero de simulacoes:\n"); scanf("%d", &sim); printf("\n"); for(j=0;j<sim;j++){ printf("Numero de iteracoes n:\n"); scanf("%d", &n); /*caso queira entrar outro valor para b printf("\n"); printf("Valor de b - qtde de dias (30 dias)\n"); scanf("%lf",&b);*/ printf("\n"); h=(b-a)/n; printf("h=%f", h); printf("\n"); printf("\n"); /*método de Euler implicito*/ for (i=0;i<=n;i++){ t[i+1]=t[i]+h; F[i]=Y(t[i]); /*metodo de Euler para calcular o ponto inicial do metodo de Newton*/ k1= f_euler(t[i],z[i]); z[i+1]=z[i]+h*k1; y[i+1]=newton( z[i], y[i], h, 0.0001, MAX); } /*fim do método*/ /*IMPRESSAO*/ printf(" i ti sol. aprox sol. exata erro global\n"); printf("\n"); //i=i-1;//usa na validação i=0; //usa na impressao da tabela while(i<=n){ //usa na impressao da tabela printf("%5d %10.3f %18.12f %18.12f %18.12f\n",i,t[i], y[i], F[i], fabs(F[i]-y[i])); i=i+1; //usa na impressao da tabela } //usa na impressao da tabela printf("\n"); val[j]=fabs(F[i-1]-y[i-1]); printf("\n"); /*validação ->|sol. aprox(h)-sol. exata|/|sol. aprox(h/2)-sol. exata|*/ //printf("Ordem de Consistencia p[%d]=%10.12f\n",j, log2(val[j-1]/val[j]));//usa na validação printf("\n"); printf("\n"); } /*usado para verificação do tempo de execução tf = clock(); tempo_gasto = ( (double) (tf - to) ) / CLOCKS_PER_SEC; printf("Tempo gasto: %lf s\n", tempo_gasto); */ printf("\n"); system("pause"); return 0; }

/*MÉTODO DE RUNGE-KUTTA-FEHLBERG*/ /* --------------------------------------------------------------------------- REFERÊNCIAS: Algorithm translated to C by: Dr. Norman Fahrer IBM and Macintosh verification by: Daniel Mathews NUMERICAL METHODS: C Programs, (c) John H. Mathews 1995

25

To accompany the text: NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992 Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A. Prentice Hall, Inc.; USA, Canada, Mexico ISBN 0-13-624990-6 Prentice Hall, International Editions: ISBN 0-13-625047-5 This free software is compliments of the author. E-mail address: in%"[email protected]" Algorithm 9.5 (Runge-Kutta-Fehlberg Method RKF45). Section 9.5, Runge-Kutta Methods, Page 461 --------------------------------------------------------------------------- */ /* --------------------------------------------------------------------------- (Metodo de Runge-Kutta-Fehlberg(RKF45)). Aproximar a solucao do problema de valor inicial y' = f(t,y) com y(a) = y_0 no intervalo [a,b] com controle do erro e tamanho variavel de passo. --------------------------------------------------------------------------- */ #include<stdio.h> #include<stdlib.h> #include<math.h> #include <time.h> #define MAX 50000 double m=pow(10,5); double k=2*pow(10,-6); double C=99*pow(10,-5); double t0=0; double y0=pow(10,3); /*Função*/ double f(double t, double y){ double Y; Y=k*(m-y)*y; return Y; } /*solucao exata*/ double Y(double t){ double y; y=m*exp(k*m*t)/(exp(k*m*t)+m*C); return y; } /* PROGRAMA PRINCIPAL */ int main() { double tol = 1E-04; /* Tolerancia do erro de discretização */ int i, j; /* contadores */ int n=100; /* número de iterações */ double a=0, b=30, y[MAX]; /* Endpoints and inital value */ double h; /* tamanho do paso */ double t[MAX]; double F[MAX]; /*solucao exata*/ double k1, k2, k3, k4, k5, k6; /* fuções a serem avaliadas */

26

double auxt, auxy; double Hmin, Hmax; /* tamanho mínimo e máximo do passo */ double Erro, S; /* Erro e constante de avaliacao do passo */ /*usado na verificação do tempo de execução clock_t to, tf; double tempo_gasto; to = clock();*/ //printf("entre com o intervalo [a,b]:\n"); //printf("por exemplo : 0, 1.4\n"); //scanf("%lf, %lf", &a, &b); //printf("intervalo escolhido: [%lf, %lf]\n", a, b); //printf("---------------------------------------------\n"); //printf("Entre com o numero de passos n\n"); //scanf("%d",&n); if(n > MAX) { printf(" programa nao preparado para mais que %d passos. FIM!\n",MAX); exit(0); } printf("-----------------------------------------\n"); /*calculando passo inicial*/ h = (b - a) / n; printf("h inicial=%f", h); printf("\n"); printf("\n"); Hmin = h / 64.0; Hmax = h * 64.0; /*condicoes iniciais*/ t[0]=t0; y[0]=y0; j= 0; while(t[j] < b) { if( (t[j] + h) > b ) h = b - t[j]; /* o ultimo passo*/ auxt = t[j]; auxy = y[j]; /* calculo das avaliacoes de f */ k1 = h * f(auxt,auxy); k2 = h * f(auxt + 0.25 * h, auxy + 0.25 * k1); k3 = h * f(auxt + 3.0*h/8.0, auxy + 3.0*k1/32.0 + 9.0*k2/32.0); k4 = h * f(auxt + 12.0*h/13.0, auxy + 1932.0*k1/2197.0 - 7200.0*k2/2197.0 + 7296.0*k3/2197.0); k5 = h * f(auxt + h, auxy + 439.0*k1/216.0 - 8.0*k2 + 3680.0*k3/513.0 - 845.0*k4/4104.0); k6 = h * f(auxt + 0.5*h, auxy - 8.0*k1/27.0 + 2.0*k2 - 3544.0*k3/2565.0 + 1859.0*k4/4104.0 - 11.0*k5/40.0); /* estimativa do erro local de truncamento */ Erro = fabs( k1/360.0 - 128.0*k3/4275.0 - 2197.0*k4/75240.0 + k5/50.0 + 2.0*k6/55.0 );

27

/* avaliacao do erro */ if ( (Erro < tol) || (h <= 2.0*Hmin) ) { F[j]=Y(t[j]); y[j+1] = auxy + 25.0*k1/216.0 + 1408.0*k3/2565.0 + 2197.0*k4/4104.0 - k5/5.0; t[j+1] = auxt + h; j++; } /* caso de divisão por zero */ if( Erro == 0 ) S = 0; /*calcula a constante S para avaliar o novo passo h*/ else S = 0.84 * sqrt( sqrt(tol * h / Erro) ); /* reducao ou aumento do passo */ if ( (S < 0.75) && (h > 2.0*Hmin) ) h /= 2.0; else if ( (S > 1.50) && (2.0 * h < Hmax) ) h *= 2.0; printf("%5d %10.3f\n",j, h); } /* fim do loop while */ /*calculo da última solucao exata*/ F[j]=Y(t[j]); printf("\n"); /* impressao da tabela de saida */ printf(" i ti sol. aprox sol. exata erro global\n"); printf("\n"); for(i=0;i<=j;i++){ printf("%5d %10.3f %18.12f %18.12f %18.12f\n",i,t[i], y[i], F[i], fabs(F[i]-y[i])); } printf("\n"); printf("\n"); /*usado para a verificação do tempo de execução tf = clock(); tempo_gasto = ( (double) (tf - to) ) / CLOCKS_PER_SEC; printf("Tempo gasto: %lf s\n", tempo_gasto); */ system("pause"); return 0; } /*(Fim do Metodo de Runge-Kutta-Fehlberg(RKF45)).*/

/*MÉTODO DE ADAMS-BASHFORTH DE 4 PASSOS*/ /*(Metodo de Adams-Bashforth de 4 Passos e Quarta Ordem). Aproximar a solucao do problema de valor inicial y' = f(t,y) com y(a) = y_0 no intervalo [a,b].

28

--------------------------------------------------------------------------- */ #include<stdio.h> #include<stdlib.h> #include<math.h> #include <time.h> #define MAX 50000 double m=pow(10,5); double k=2*pow(10,-6); double C=99*pow(10,-5); double t0=0; double y0=pow(10,3); /*Função*/ double f(double t, double y){ double Y; Y=k*(m-y)*y; return Y; } /*solucao exata*/ double Y(double t){ double y; y=m*exp(k*m*t)/(exp(k*m*t)+m*C); return y; } /* PROGRAMA PRINCIPAL */ int main() { int i, n, j, sim; /*usado na verificação do tempo de execução clock_t to, tf; double tempo_gasto; to = clock();*/ double t[MAX], y[MAX]/*sol. aprox.*/, h, F[MAX]/*sol. exata*/, val[MAX], a, b=30, k1,k2,k3,k4; a=0; t[0]=t0;/*condição inicial*/ y[0]=y0;/*condição inicial*/ printf("Numero de simulacoes:\n"); scanf("%d", &sim); printf("\n"); for(j=0;j<sim;j++){ printf("Numero de iteracoes n:\n"); scanf("%d", &n); printf("\n"); //printf("Valor de b - qtde de dias (30 dias)\n"); //scanf("%lf",&b); h=(b-a)/n; printf("h=%f", h); printf("\n"); printf("\n"); /*método RK44 (classico)- calculo de y[1], y[2] e y[3]*/ for (i=0;i<=3;i++){ t[i+1]=t[i]+h; F[i]=Y(t[i]); k1= f(t[i],y[i]); k2= f(t[i]+0.5*h,y[i]+0.5*h*k1);

29

k3= f(t[i]+0.5*h,y[i]+0.5*h*k2); k4= f(t[i]+h,y[i]+h*k3); y[i+1]=y[i]+h*(k1+2*k2+2*k3+k4)/6.0; } /* Adams-bashforth */ for (i=3;i<=n;i++){ t[i+1]=t[i]+h; F[i]=Y(t[i]); y[i+1]=y[i]+h*(55*f(t[i],y[i])-59*f(t[i-1],y[i-1])+37*f(t[i-2],y[i-2])-9*f(t[i-3],y[i-3]))/24; } /*fim do método*/ //i=i-1;// usado para fazer a validacao printf(" i ti sol. aprox sol. exata erro global\n"); printf("\n"); for(i=0;i<=n;i++){ //for utilizado para a impressao da tabela /*impressão da tabela*/ printf("%5d %10.3f %18.12f %18.12f %18.12f\n",i,t[i], y[i], F[i], fabs(F[i]-y[i])); } printf("\n"); val[j]=fabs(F[i-1]-y[i-1]); printf("\n"); /*validação ->|sol. aprox(h)-sol. exata|/|sol. aprox(h/2)-sol. exata|*/ //printf("Ordem de Consistencia p[%d]=%10.12f\n",j, log2(val[j-1]/val[j-0])); //printf("\n"); } /*usado na verificação do tempo de execução tf = clock(); tempo_gasto = ( (double) (tf - to) ) / CLOCKS_PER_SEC; printf("Tempo gasto: %lf s\n", tempo_gasto);*/ printf("\n"); system("pause"); return 0; }

/*MÉTODO DE ADAMS-MOULTON DE 4 PASSOS*/ /*(Metodo de Adams-Moulton de 4 Passos e Quinta Ordem). Aproximar a solucao do problema de valor inicial y' = f(t,y) com y(a) = y_0 no intervalo [a,b]. --------------------------------------------------------------------------- */ #include<stdio.h> #include<stdlib.h> #include<math.h> #include <time.h> #define MAX 50000 double m=pow(10,5); double k=2*pow(10,-6); double C=99*pow(10,-5); double t0=0; double y0=pow(10,3);

30

/*Função*/ double f(double t, double y){ double Y; Y=k*(m-y)*y; return Y; } /*Função a ser avaliada no EULER para o cálculo do ponto inicial no metodo de Newton*/ double f_euler(double tz, double z){ double Y; Y=k*(m-z)*z; return Y; } /*solução exata*/ double Y(double t){ double y; y=m*exp(k*m*t)/(exp(k*m*t)+m*C); return y; } /*FUNCAO IMPLÍCITA*/ double f1(double x0, double yk, double yk1, double yk2, double yk3, double h){ double x; x=-x0+yk+h*251*k*x0*(m-x0)/720+h*(646*k*yk*(m-yk) - 264*k*yk1*(m-yk1) + 106*k*yk2*(m-yk2) - 19*k*yk3*(m-yk3))/720; return (x); } /*DERIVADA DA FUNCAO IMPLICITA*/ double Df1(double x0, double h){ double dx; dx=-h*251*x0/360 + h*251*k*m/720 - 1; return (dx); } /*METODO DE NEWTON*/ double newton(double x0, double yk, double yk1, double yk2, double yk3 , double h, double tol, int N){ int i; double x; for(i=0;1<=N;i++){ x=x0-f1(x0, yk, yk1, yk2, yk3, h)/Df1(x0,h); if(fabs(x-x0)<tol) return (x); else x0=x; } printf("o metodo falhou apos N="); return (N); } /* PROGRAMA PRINCIPAL */ int main() { int i, n, j, sim; /*usado para calcular o tempo de execução do programa clock_t to, tf; double tempo_gasto; to = clock();*/ double t[MAX], y[MAX]/*sol. aprox.*/, h, F[MAX]/*sol. exata*/, val[MAX], a, b=30, z[MAX], w1, k1,k2,k3,k4,k5,k6 ; a=0; t[0]=t0;/*condição inicial*/ y[0]=y0;/*condição inicial*/ printf("Numero de simulacoes:\n"); scanf("%d", &sim);

31

printf("\n"); for(j=0;j<sim;j++){ printf("Numero de iteracoes n:\n"); scanf("%d", &n); printf("\n"); //printf("Valor de b - qtde de dias (30 dias)\n"); //scanf("%lf",&b); h=(b-a)/n; printf("h=%f", h); printf("\n"); printf("\n"); /*método RK56 - calculo de y[1], y[2] e y[3]*/ for (i=0;i<=3;i++){ t[i+1]=t[i]+h; F[i]=Y(t[i]); k1 = f(t[i] , y[i]); k2 = f(t[i] + h/4 , y[i] + h*k1/4); k3 = f(t[i] + 3*h/8, y[i] + 3*h*k1/32 + 9*h*k2/32); k4 = f(t[i] + 12*h/13 , y[i] + 1932*h*k1/2197 - 7200*h*k2/2197 + 7296*h*k3/2197 ); k5 = f(t[i] + h , y[i] + 439*h*k1/216 - 8*h*k2 + 3680*h*k3/513 - 845*h*k4/4104); k6 = f(t[i] + h/2 , y[i] - 8*h*k1/27 + 2*h*k2 - 3544*h*k3/2565 + 1859*h*k4/4104 - 11*h*k5/40); y[i+1]=y[i]+h*(16*k1/135 + 6656*k3/12825 + 28561*k4/56430 - 9*k5/50 + 2*k6/55); } /* Adams-Moulton */ for (i=3;i<=n;i++){ t[i+1]=t[i]+h; F[i]=Y(t[i]); /*metodo de Euler para calcular o ponto inicial do metodo de Newton*/ w1= f_euler(t[i],z[i]); z[i+1]=z[i]+h*w1; y[i+1]=newton( z[i], y[i], y[i-1], y[i-2] , y[i-3] , h, 0.00000000001, MAX); } //i=i-1;//usa na validação /*fim do método*/ printf(" i ti sol. aprox sol. exata erro global\n"); printf("\n"); for(i=0;i<=n;i++){ /*impressão da tabela*/ printf("%5d %10.3f %18.12f %18.12f %18.12f\n",i,t[i], y[i], F[i], fabs(F[i]-y[i])); } printf("\n"); printf("\n"); val[j]=fabs(F[i-1]-y[i-1]); printf("\n"); /*validação ->|sol. aprox(h)-sol. exata|/|sol. aprox(h/2)-sol. exata|*/ // printf("Ordem de Consistencia p[%d]=%10.12f\n",j, log2(val[j-1]/val[j-0])); printf("\n"); } printf("\n"); /*usado para calcular o tempo de execução tf = clock(); tempo_gasto = ( (double) (tf - to) ) / CLOCKS_PER_SEC; printf("Tempo gasto: %lf s\n", tempo_gasto); */ system("pause"); return 0; }

Documents

EXERCÍCIO COMPUTACIONAL – MAP 5725 RELATÓRIO · mÉtodo de runge-kutta-fehlberg (rkf) O método RKF utiliza os métodos RK45 e RK56, a idéia desse método é obter uma estimativa