Erro: Leapfrog
Ir para navegação
Ir para pesquisar
#include<stdio.h>
#include<math.h>
#include<string.h>
#include<stdlib.h>
void atualizar_onda(double u_new[],double u_now[],double u_old[],int jmax,double k);
void main()
{
FILE *arq;
arq = fopen("erro_leapfrog.dat", "w");
int i, j, jmax;
//tamanho da corda: jmax-1
jmax = 50;
//u: posicao da corda
/* u_old em t-1 */
/* u_now em t */
/* u_new em t+1 */
//k = dt/dx (aqui, dx = 1, k=dt)
double u_new[jmax], t, u_old[jmax], u_now[jmax], k, erro, tmax;
//tmax: tempo final
tmax = 100;
//variamos o dt para fazer um grafico do erro em função de dt
int n=1;
while( n < 1200)
{
k = 0.001*n;
//condicao inicial
for (j = 0 ; j < jmax ; j++)
{
u_now[j] = sin(M_PI*j/(jmax - 1));
}
//condicao de contorno
u_new[0] = 0;
u_new[jmax-1] = 0;
u_old[0] = 0;
u_old[jmax-1] = 0;
//calculo de u para t=-dt, utilizando o metodo de leapfrog - aqui usamos que du/dt = 0 em t=0
for (j = 1; j< jmax-1 ; j++)
{
u_old[j] = u_now[j] + 0.5 * pow(k,2) * (u_now[j+1] - 2 * u_now[j] + u_now[j-1]);
}
//atualizacao da onda
for(t = 0 ; t < tmax ; t+=k)
{
//para que o ultimo tempo em que se calcula a posicao seja tmax, independentemente do valor de dt
if(t+k >tmax)
{
k = tmax - t;
}
atualizar_onda(u_new,u_now,u_old,jmax,k);
memcpy(u_old,u_now, sizeof(double)*jmax);
memcpy(u_now,u_new, sizeof(double)*jmax);
}
//calculo do erro
erro = 0;
for(j = 0 ; j < jmax ; j++)
{
double analitica = cos((M_PI*t)/(jmax - 1))*sin((M_PI*j)/(jmax - 1));
erro += fabs(analitica - u_now[j]);
}
erro = erro/jmax;
fprintf(arq,"%.12lf %.12lf\n", n*0.001,erro);
n++;
}
fclose(arq);
}
void atualizar_onda(double u_new_[],double u_now_[],double u_old_[],int jmax,double k)
{
int j;
for(j = 1 ; j < jmax-1 ; j++)
{
u_new_[j] = 2 * u_now_[j] - u_old_[j] + pow(k,2) *( u_now_[j+1] - 2 * u_now_[j] + u_now_[j-1] );
}
}