Problem 3 main.c

#include <stdio.h>
#include <math.h>

double integrand(double x)
{
    if (x!=0){
    return x / sqrt(cosh(4) - cosh(4-x*x));
    }
    else return 1/sqrt(sinh(4));
}

double f (double x)
{
    return -sinh(x);
}

double E (double x, double v)
{
    return v*v/2 + cosh(x);
}

void derivs(double t, double y[], double dydt[])
{
    dydt[0] = y[1];
    dydt[1] = - sinh(y[0]);

}


int main()
{
double simp(), f(), a, b, T, EPS, ans, prevans;
int n=0;
a=0;
b=2;
T=0;
prevans = 1e50;
EPS = 1e-5;


for(n=2; fabs(T - prevans)>EPS; n*=2)
 {
    ans = simp(integrand,a,b,n);
    prevans = T;
    T = 4*sqrt(2)* ans;
}
   printf("T = %10.5f \n\n", T);

double f(), x, v, h, t, E(), DeltaEmax, E_exact = cosh(4), k;
h = T/50;
t = 0;
DeltaEmax=0;
printf("Verlet Algorithm:\n\n"
       " n         DeltaE_max\n\n");
for(n=1; n<=100; n*=10){
v = 0; x = 4;
    for(t=h; t<n*(T)+h; t+=h)
    {
        v+= h / 2 * f(x);
        x+= h*v;
        v+=h/2*f(x);
        k = fabs(E(x,v)-E_exact);
        if(k>DeltaEmax)
            DeltaEmax = k;

    }
printf("%3d %15.5f \n", n,  DeltaEmax);
}

printf("\n RK2 Algorithm:\n\n"
       " n          DeltaE_max\n\n");
int M = 2;
double rk2(), y[M];
extern void derivs();
DeltaEmax= 0;
for(n=1; n<=100; n*=10){
    y[0]=4; y[1] = 0;
    for(t=h; t<=n*(T)+h; t+=h){
        rk2(t, y, derivs, h, M);
        k = fabs(E(y[0],y[1])-E_exact);
        if(k>DeltaEmax)
            DeltaEmax = k;

    }
 printf("%3d %15.5f \n", n, DeltaEmax);

}




    return 0;
}