Lax-Friedrichs

#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("lax.txt", "w+");

    int tmax, 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: 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_old[0]=0;
    u_old[jmax-1]=0;

    u_new[0] = 0;
    u_new[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);
    }

	//onda no tempo t=tmax
    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;

    //para j=1 e j=jmax-2 (penultimas posicoes), necessitamos de u em x=-1 e x=jmax, respectivamente, sendo que essas posicoes ficam fora do vetor (vai de 0 a jmax-1). Aqui, foi feito u(x=-1)=-u(x=1) e u(x=jmax)=-u(x=jmax-2)
    j=1;
    u_new_[j] = u_now_[j] + 0.5 * (u_now_[j-1] + u_now_[j+1]) - 0.5 * (u_old_[j-1] + u_old_[j+1]) + 0.25 * pow(k,2) * (u_old_[j+2] - 2 * u_old_[j] - u_old_[j]);

    for(j = 2 ; j < jmax-2 ; j++)
    {
        u_new_[j] = u_now_[j] + 0.5 * (u_now_[j-1] + u_now_[j+1]) - 0.5 * (u_old_[j-1] + u_old_[j+1]) + 0.25 * pow(k,2) * (u_old_[j+2] - 2 * u_old_[j] + u_old_[j-2]);
    }

    j=jmax-2;
    u_new_[j] = u_now_[j] + 0.5 * (u_now_[j-1] + u_now_[j+1]) - 0.5 * (u_old_[j-1] + u_old_[j+1]) + 0.25 * pow(k,2) * (-u_old_[j] - 2 * u_old_[j] + u_old_[j-2]);



}