Solution Lecture b3

! Solve the heat equation with periodic boundary condition
implicit none
integer, parameter :: M = 10
real, parameter :: pi = 4. * atan(1.)
integer n, k, nmax
real h, tau, x
real, dimension(0:M-1) :: u, v

h = 1. / M
tau = h * h / 4.
nmax = 0.125 / tau

do k=0, M-1
    x = k * h
    u(k) = cos(2 * pi * x)
    print *, 0., x, 0.
end do
print *

do n=0, nmax-1
    do k = 0, M-1
        v(k) = u(k) + (tau / (h * h)) * (u(mod(k + M - 1, M)) &
        - 2 * u(k) + u(mod(k + 1, M)))
    end do

    do k = 0, M-1
        print *, (n + 1) * tau, k * h, abs(v(k) - exp(-4*pi**2*(n+1)*tau) * cos(2. * pi * k * h))
    end do
    print *

    u = v
end do
end

implicit none
integer, parameter :: M = 10
integer n, k, nmax
real h, tau, x, a, b
real, dimension(0:M) :: u, v, error

h = 1. / M
tau = h * h / 4.
nmax = 0.5 / tau

! Initial data u_0 = 0
do k = 0, M
    x = k * h
    u(k) = x * (1 - x)
enddo
! Print the initial error (0)
 call printstep(0., 0. * u)

a = 1
b = -1

do n=0, nmax-1
    v(0) = u(0) + 2 * tau / (h * h) * (u(1) - u(0) - h * a)
    do k = 1, M-1
        v(k) = u(k) + (tau / (h * h)) * (u(k - 1) - 2 * u(k) + u(k + 1))
    end do
    v(M) = u(M) + 2 * tau / (h * h) * (u(M-1) - u(M) + h * b)

    do k = 0, M
        x = k * h
        error(k) = abs(v(k) - x * (1 - x) + 2 * (n + 1) * tau)
    enddo
    call printstep((n+1) * tau, error)

    u = v
end do
 contains
    ! Prints the values of the solution at a given time step
    subroutine printstep(tn, un)
        real tn
        real, dimension(0:M) :: un
        integer k
        do k = 0, M
            print *, tn, k * h, un(k)
        end do
        ! print an empty line (needed for gnuplot)
        print *
    end subroutine
end

implicit none
integer, parameter :: M = 10
integer n, k, nmax
real h, tau, x, f
real, dimension(0:M) :: u, v

h = 1. / M
tau = h * h / 4.
nmax = 0.5 / tau

! Initial data u_0 = 0
u = 0.
! Print the initial data
 call printstep(0., u)

do n=0, nmax-1
    v(0) = u(0)
    do k = 1, M-1
        f = 10 - 20 * abs(k * h - 0.5)
        v(k) = u(k) + (tau / (h * h)) * (u(k - 1) - 2 * u(k) + u(k + 1)) + tau * f
    end do
    v(M) = u(M)

    call printstep((n+1) * tau, v)

    u = v
end do
 contains
    ! Prints the values of the solution at a given time step
    subroutine printstep(tn, un)
        real tn
        real, dimension(0:M) :: un
        integer k
        do k = 0, M
            print *, tn, k * h, un(k)
        end do
        ! print an empty line (needed for gnuplot)
        print *
    end subroutine
end