Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #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;
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement