Leapfrog

#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("leapfrog.txt", "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
    double u_new[jmax], t, u_old[jmax], u_now[jmax], k, erro, tmax;

    //tmax: tempo final
    tmax = 100;
    k = 0.2;

    //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)
    {


        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);
    }

    for(j=0; j<jmax; j++)
    {
        fprintf(arq,"%d %lf\n", j, u_new[j]);
    }


    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] );
    }


}