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
- Write a do loop that sums the squares of the first
nintegers (from 1 ton). - Define a new function
sum_of_squaresthat calculates this sum and returns it. - In the
mainprogram, callsum_of_squaresand 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
- Create a function named factorial that computes the factorial of a non-negative integer using a do loop.
- 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:
- Given the time , compute the new position at time :
- Compute the new acceleration:
- Compute the new velocity:
Here a to do list that can help:
-
Define Parameters: define the mass
m, spring constantk, time stepdt, and total simulationtime. -
Set Up Initial Conditions: det the initial position
x_0and initial velocityv_0. -
Implement the Verlet Integration Algorithm: use the following equations to update the position and velocity at each time step:
where .
-
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()andrand():srand(time(0)); // Seed the random number generator int r = rand() % 100; // Generates a random number between 0 and 99where:
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 andRAND_MAX(whereRAND_MAXis 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()andrandom_number():real :: r call random_seed() ! Seed the random number generator call random_number(r) ! Generates a random number between 0 and 1where:
call random_seed(): Initializes the random number generator (typically based on the system clock if no seed is provided).call random_number(x): Fillsxwith 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:
- 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:
- Calculate the total energy: for each particle, calculate the total energy based on the generated velocities:
where is the magnitude of the momentum .
-
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 .
-
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
-
https://en.wikipedia.org/wiki/Verlet_integration ↩