Exercises related to C_BASIC

Here some exercises with the solutions in C and FORTRAN. Feel free to start from them and expand them as you wish. For example, for some of them can be usefull to print the output!

Exercise 1: Sum of the first squares

  1. Write a do loop that sums the squares of the first n integers (from 1 to n).
  2. Define a new function sum_of_squares that calculates this sum and returns it.
  3. In the main program, call sum_of_squares and print the result.

Solution

C code
#include <stdio.h>

// Function to calculate the sum of squares from 1 to n
int sum_of_squares(int n) {
    int sum = 0;
    for (int i = 1; i <= n; i++) {
        sum += i * i;  // Add the square of i to sum
    }
    return sum;
}

int main() {
    int n;

    // Prompt the user for input
    printf("Enter a positive integer n: ");
    scanf("%d", &n);

    // Call the function and store the result
    int result = sum_of_squares(n);

    // Output the result
    printf("The sum of the squares of the first %d integers is: %d\n", n, result);
    return 0;
}
FORTRAN code
program sum_of_squares_program
    implicit none

    integer :: n, result

    ! Prompt the user for input
    print *, "Enter a positive integer n: "
    read(*, *) n

    ! Call the function and store the result
    result = sum_of_squares(n)

    ! Output the result
    print *, "The sum of the squares of the first", n, "integers is:", result

contains

    ! Function to calculate the sum of squares from 1 to n
    integer function sum_of_squares(n)
        implicit none
        integer :: n, i
        integer :: sum

        sum = 0
        do i = 1, n
            sum = sum + i * i  ! Add the square of i to sum
        end do

        sum_of_squares = sum  ! Return the computed sum
    end function sum_of_squares

end program sum_of_squares_program

Exercise 2: Factorial

  1. Create a function named factorial that computes the factorial of a non-negative integer using a do loop.
  2. In the main program, use the factorial function to compute the factorial of numbers from 1 to 10 and print the results.

Solution

C code
#include <stdio.h>

// Function to compute factorial
int factorial(int n) {
    int result = 1;
    for (int i = 1; i <= n; i++) {
        result *= i;  // Multiply result by i
    }
    return result;
}

int main() {
    // Print factorials of numbers from 1 to 10
    printf("Factorials from 1 to 10:\n");
    for (int i = 1; i <= 10; i++) {
        printf("%d! = %llu\n", i, factorial(i)); // Call factorial function
    }
    return 0;
}
FORTRAN code
program factorial_program
    implicit none

    integer :: i, n
    print *, "Factorials from 1 to 10:"

    ! Loop through numbers from 1 to 10
    do i = 1, 10
        print *, i, "! =", factorial(i)
    end do

contains
    ! Function to compute factorial
    integer function factorial(n)
        implicit none
        integer :: n, result, j

        result = 1
        do j = 1, n
            result = result * j   ! Multiply result by j
        end do

        factorial = result  ! Return the computed factorial
    end function factorial
end program factorial_program

Exercise 3: Mechanical Energy

Write a program that simulates the motion of a point mass ( kg) under the influence of gravity. The program should calculate and print the height (), velocity (), kinetic energy (), and potential energy () of the mass as functions of time. Assume that the object starts from an initial height m and with an initial velocity m/s. Stop the code when the point touches the ground.

Formulae

  • Height as a function of time:
  • Velocity as a function of time:
  • Kinetic Energy:
  • Potential Energy:

Solution

C code
#include <stdio.h>

#define g 9.81  // Gravitational acceleration (m/s^2)

// Function to compute kinetic energy
double kinetic_energy(double mass, double velocity) {
    return 0.5 * mass * velocity * velocity;
}

// Function to compute potential energy
double potential_energy(double mass, double height) {
    return mass * g * height;
}

// Function to compute velocity at time t
double velocity_at_time(double v0, double t) {
    return v0 - g * t;  // v(t) = v0 - g * t
}

// Function to compute height at time t
double height_at_time(double h0, double v0, double t) {
    return h0 + v0 * t - 0.5 * g * t * t;  // h(t) = h0 + v0 * t - 1/2 * g * t^2
}

int main() {
    double t = 0.0;            // Time starts at 0
    double dt = 0.1;           // Time step
    double total_time = 10.0;  // Total simulation time
    double m  = 1.0;           // Mass of the point (kg)
    double v0 = 10;            // Initial velocity (m/s) 
    double h0 = 100;           // Initial height (m)

    // Header for the output
    printf("\nTime(s)\t\tKinetic Energy(J)\tPotential Energy(J)\n");
    printf("------------------------------------------------------------\n");

    // Loop over time and compute KE and PE
    for (t=0; t <= total_time; t+=dt) {
        double h = height_at_time(h0, v0, t);  // Height at time t
        double v = velocity_at_time(v0, t);    // Velocity at time t

        if (h < 0) h = 0;  // Stop the object from going below ground level

        double KE = kinetic_energy(m, v);      // Compute kinetic energy
        double PE = potential_energy(m, h);    // Compute potential energy

        // Print time, kinetic energy, and potential energy
        printf("%.2lf\t\t%.2lf\t\t\t%.2lf\n", t, KE, PE);

        // Stop the loop if the object hits the ground (h = 0)
        if (h == 0 && v <= 0) break;
    }

    return 0;
}
FORTRAN code
program energy_simulation
  implicit none

  ! Constants
  real, parameter :: g = 9.81   ! Gravitational acceleration (m/s^2)

  ! Variables
  real :: t, dt, total_time, m, v0, h0, h, v, KE, PE

  ! Initialize values
  t = 0.0                     ! Time starts at 0
  dt = 0.1                     ! Time step
  total_time = 10.0            ! Total simulation time
  m = 1.0                      ! Mass of the point (kg)
  v0 = 10.0                    ! Initial velocity (m/s)
  h0 = 100.0                   ! Initial height (m)

  ! Output header
  print *, 'Time(s)', '    Kinetic Energy(J)', '   Potential Energy(J)'
  print *, '---------------------------------------------------------'

  ! Loop over time and compute KE and PE
  do while (t <= total_time)
    
    ! Compute height and velocity at time t
    h = height_at_time(h0, v0, t)
    v = velocity_at_time(v0, t)

    ! Prevent object from going below ground level
    if (h < 0.0) h = 0.0

    ! Compute kinetic and potential energy
    KE = kinetic_energy(m, v)
    PE = potential_energy(m, h)

    ! Print the results
    print "(F6.2, 2X, F10.2, 3X, F10.2)", t, KE, PE

    ! Stop the loop if the object hits the ground and velocity is downward
    if (h == 0.0 .and. v <= 0.0) exit

    ! Increment time
    t = t + dt

  end do

contains

  ! Function to compute kinetic energy
  real function kinetic_energy(mass, velocity)
    real, intent(in) :: mass, velocity
    kinetic_energy = 0.5 * mass * velocity * velocity
  end function kinetic_energy

  ! Function to compute potential energy
  real function potential_energy(mass, height)
    real, intent(in) :: mass, height
    potential_energy = mass * g * height
  end function potential_energy

  ! Function to compute velocity at time t
  real function velocity_at_time(v0, t)
    real, intent(in) :: v0, t
    velocity_at_time = v0 - g * t
  end function velocity_at_time

  ! Function to compute height at time t
  real function height_at_time(h0, v0, t)
    real, intent(in) :: h0, v0, t
    height_at_time = h0 + v0 * t - 0.5 * g * t * t
  end function height_at_time

end program energy_simulation

Exercise 4: Simulation of a Simple Harmonic Oscillator

A simple harmonic oscillator is described by the equation of motion:

where:

  • is the mass of the oscillator,
  • is the spring constant,
  • is the displacement from the equilibrium position.

Your task is to numerically solve this differential equation using the velocity Verlet integration method, which is particularly useful for simulating physical systems.

Velocity Verlet 1

For a second-order differential equation of the type

with initial conditions

an approximate numerical solution at the time with time step can be obtained by the following method:

  1. Given the time , compute the new position at time :
  2. Compute the new acceleration:
  3. Compute the new velocity:

Here a to do list that can help:

  1. Define Parameters: define the mass m , spring constant k , time step dt , and total simulation time.

  2. Set Up Initial Conditions: det the initial position x_0 and initial velocity v_0 .

  3. Implement the Verlet Integration Algorithm: use the following equations to update the position and velocity at each time step:

    where .

  4. Output the Results: print the position and velocity of the oscillator at each time step.

Solution

C code
#include <stdio.h>

void verlet_integration(double m, double k, double x0, double v0, double dt, double total_time) {
    double x = x0; // current position
    double v = v0; // current velocity
    double a = -k/m * x; // initial acceleration
    double time = 0.0;

    printf("Time\tPosition\tVelocity\n");

    while (time <= total_time) {
        printf("%.2f\t%.4f\t%.4f\n", time, x, v);

        double x_new = x + v * dt + 0.5 * a * dt * dt; // update position
        double a_new = -k/m * x_new; // calculate new acceleration
        v = v + 0.5 * (a + a_new) * dt; // update velocity
        x = x_new; // update position for next iteration
        a = a_new; // update acceleration

        time += dt; // increment time
    }
}

int main() {
    double m = 1.0;      // mass in kg
    double k = 10.0;     // spring constant in N/m
    double dt = 0.01;    // time step in seconds
    double total_time = 2.0; // total simulation time in seconds
    double x0 = 1.0;     // initial position in meters
    double v0 = 0.0;     // initial velocity in m/s

    verlet_integration(m, k, x0, v0, dt, total_time);

    return 0;
}
FORTRAN code
program verlet_integration
    implicit none

    real(8) :: m, k, x, v, a, dt, total_time, time, x_new, a_new, x0, v0

    ! Define parameters
    m = 1.0d0           ! Mass in kg
    k = 10.0d0          ! Spring constant in N/m
    dt = 0.01d0         ! Time step in seconds
    total_time = 2.0d0  ! Total simulation time in seconds
    x0 = 1.0d0          ! Initial position in meters
    v0 = 0.0d0          ! Initial velocity in m/s

    ! Initial conditions
    x = x0              ! Current position
    v = v0              ! Current velocity
    a = -k/m * x        ! Initial acceleration
    time = 0.0d0        ! Start time

    ! Print header
    print *, "Time", "Position", "Velocity"

    ! Time loop using Verlet integration
    do while (time <= total_time)
        ! Print current state
        print '(F6.2, 2X, F8.4, 2X, F8.4)', time, x, v

        ! Compute new position
        x_new = x + v * dt + 0.5d0 * a * dt * dt

        ! Compute new acceleration based on new position
        a_new = -k/m * x_new

        ! Update velocity
        v = v + 0.5d0 * (a + a_new) * dt

        ! Update position and acceleration for the next step
        x = x_new
        a = a_new

        ! Increment time
        time = time + dt
    end do

end program verlet_integration

Exercise 5: Random Numbers

Write a program that generates N random integers between 1 and 100, where N is specified by the user. The program should print each number as it is generated, and at the end, it should calculate and display the sum and average of the numbers. Use srand() to ensure different random numbers are generated each time the program runs.

[!TIP]

C

To gerenate random numbers unfiromly distributed in a given range, we can use srand() and rand():

srand(time(0));       // Seed the random number generator
int r = rand() % 100; // Generates a random number between 0 and 99

where:

  • srand(unsigned int seed): Initializes the random number generator. Using srand(time(0)) seeds it with the current time, ensuring different random numbers on each run.
  • rand(): Generates a pseudo-random integer between 0 and RAND_MAX (where RAND_MAX is the macro that is defined in the <stdlib.h> library). To limit the range, use modulo: rand() % max + min.

FORTRAN

To gerenate random numbers unfiromly distributed in a given range, we can use random_seed() and random_number():

real :: r
call random_seed()      ! Seed the random number generator
call random_number(r)   ! Generates a random number between 0 and 1

where:

  • call random_seed(): Initializes the random number generator (typically based on the system clock if no seed is provided).
  • call random_number(x): Fills x with a random number (or array of random numbers) uniformly distributed between 0 and 1.

Solution

C code
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

int main() {
    int N, i, num, sum = 0;
    float average;

    // Ask the user for the number of random numbers (N)
    printf("Enter the number of random numbers to generate: ");
    scanf("%d", &N);

    // Seed the random number generator with the current time
    srand(time(0));

    printf("\nRandom numbers between 1 and 100:\n");

    // Generate N random numbers between 1 and 100
    for (i = 0; i < N; i++) {
        num = (rand() % 100) + 1;
        printf("%d ", num);
        sum += num;  // Add the number to the sum
    }

    // Calculate the average
    average = (float)sum / N;

    printf("\n\nSum of the numbers: %d\n", sum);
    printf("Average of the numbers: %.2f\n", average);

    return 0;
}
FORTRAN code
program random_numbers
  implicit none

  integer :: N, i, num, sum
  real :: average
  real :: r

  ! Ask the user for the number of random numbers (N)
  print *, "Enter the number of random numbers to generate: "
  read *, N

  ! Seed the random number generator with the current time
  call random_seed()

  sum = 0
  print *, "Random numbers between 1 and 100:"

  ! Generate N random numbers between 1 and 100
  do i = 1, N
      call random_number(r)            ! Generates a random number between 0 and 1
      num = int(r * 100) + 1           ! Scale the number to the range [1, 100]
      print *, num
      sum = sum + num                  ! Add the number to the sum
  end do

  ! Calculate the average
  average = real(sum) / real(N)

  print *, " "
  print *, "Sum of the numbers: ", sum
  print *, "Average of the numbers: ", average
end program random_numbers

Exercise 6: Maxwell-Boltzmann Distribution and Energy Calculation

Consider particles with mass kg in a three-dimensional space. The Hamiltonian of the single particle is

where is the momentum and and it is distributed according to the Maxwell-Boltzmann distribution:

where J/K is the Boltzmann constant.

The task is to compute the energy of the system and to compare the result with the equipartition theorem: .

Tasks:

  1. Generate the momenta: For a system of particles at a temperature K, generate the velocity components , , from a normal distribution with a standard deviation given by:

  1. Calculate the total energy: for each particle, calculate the total energy based on the generated velocities:

where is the magnitude of the momentum .

  1. Compare to the equipartition theorem: compare the calculated energy to the value predicted by the equipartition theorem, which states that the total energy of the system should be .

  2. Vary the number of particles: Perform the simulation for different values of (e.g., 10, 100, 1000, 10000, etc.) and observe how the results converge to the theoretical value as N increases.

Solution

C code
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

// Constants
const double k_B = 1.38e-23;      // Boltzmann constant (J/K)
const double T = 300.0;           // Temperature in Kelvin
const double mass = 1.e-27;       // Mass of the particle (arbitrary units)

// Function to generate normally distributed random numbers
double random_normal(double mean, double stddev) {
    double u1 = (double)rand() / RAND_MAX;
    double u2 = (double)rand() / RAND_MAX;
    return mean + stddev * sqrt(-2.0 * log(u1)) * cos(2.0 * M_PI * u2);
}

int main() {
    int N;            // Number of particles
    double vx, vy, vz, p, energy, sigma;
    int i, j;

    // Seed for random number generation
    srand(time(NULL));

    // Standard deviation for the normal distribution
    sigma = sqrt(k_B * T / mass);

    N = 1;
    for (i = 1; i <= 5; i++) {
        N *= 10; // Increase the number of particles each loop
        energy = 0.0;

        // Generate momentum values according to Maxwell-Boltzmann distribution
        for (j = 0; j < N; j++) {
            vx = random_normal(0.0, sigma);
            vy = random_normal(0.0, sigma);
            vz = random_normal(0.0, sigma);
            p = mass * sqrt(vx * vx + vy * vy + vz * vz);
            energy += (p * p) / (2.0 * mass);
        }

        // Print the results
        printf("----------------------------------------------------\n");
        printf("Number of particles: %d\n", N);
        printf("Energy of the system: %.5e J\n", energy);
        printf("Equipartition theorem (U=3/2NKT): %.5e J\n", (3.0 / 2.0) * N * k_B * T);
        printf("Relative difference: %.5f%%\n", (fabs(energy - (3.0 / 2.0) * N * k_B * T) * 2.0 / (energy + (3.0 / 2.0) * N * k_B * T)) * 100.0);
    }

    return 0;
}
FORTRAN code
program average_energy_mb
  implicit none

  integer            :: N             ! Number of particles
  real(8), parameter :: mass = 1.e-27 ! Mass of the particle (arbitrary units)
  real(8), parameter :: k = 1.38e-23  ! Boltzmann constant (arbitrary units)
  real(8), parameter :: T = 300.0     ! Temperature in Kelvin
  real(8)            :: p, energy, sigma, vx,vy,vz
  integer            :: i,j

  ! Seed for random number generation
  call random_seed()

  ! Standard deviation for the normal distribution
  sigma = sqrt(k * T / mass)

  N = 1
  do i = 1, 5 !Change the number of particles
    N = N * 10

    ! Generate momentum values according to Maxwell-Boltzmann distribution
    energy = 0.
    do j = 1, N
      vx = random_normal(0.0d0, sigma)
      vy = random_normal(0.0d0, sigma)
      vz = random_normal(0.0d0, sigma)
      p  = mass*sqrt(vx*vx+vy*vy+vz*vz)
      energy = energy + (p**2/ (2.0d0 * mass))
    end do

    ! Print the results
    print *, "----------------------------------------------------"
    print '(A35,I12)', "Number of particles:", N
    print '(A35,E12.5,A)', "Energy of the system:", energy, " J"
    print '(A35,E12.5,A)', "Equipartition theorem (U=3/2NKT):", 3./2.*N*K*T, " J"
    print '(A35,E12.5,A)', "Relative difference:",(abs(energy-(3./2.*N*K*T))*2./(energy+(3./2.*N*K*T)))*100., "%"
  enddo

contains

  ! Subroutine to generate normally distributed random numbers
  function random_normal(mean, stddev)
      double precision :: random_normal, mean, stddev
      double precision :: u1, u2
      u1 = rand()
      u2 = rand()
      random_normal = mean + stddev * sqrt(-2.0d0 * log(u1)) * cos(2.0d0 * 3.14159265358979d0 * u2)
  end function

end program

References


  1. https://en.wikipedia.org/wiki/Verlet_integration