22) HW2 Review#
Last time#
Solving Systems
Direct methods
Iterative methods
Example of a PDE
Matrix norm
Matrix condition number
Today#
Review of HW2
Submission expectations
1. Review of HW2#
Parts 1-3#
Here is my solution for the bisection.f90 code:
! bisection.f90 program
program main
use params
implicit none
interface
subroutine bisect_hist(hist, a, b, tol, it)
real, dimension(:), intent(inout) :: hist
real, intent(inout) :: a, b
real, intent(in) :: tol
integer, intent(out) :: it
end subroutine bisect_hist
end interface
real :: tol
integer :: niter1, niter2
real :: r1, r2
real, external :: f
INTEGER :: i
real :: hist1(n), hist2(n)
tol = 1.0e-4
call bisect_hist(hist1, a1, b1, tol, niter1)
if (niter1 >= 1) then
write(*,'(A,1X,F12.7,1X,A,1X,ES12.4)') "Root on [-4,-2]:", hist1(niter1-1), "f(x*)=", f(hist1(niter1-1))
write(*,'(A,1X,I0)') "Iterations:", niter1
OPEN(1,file='output.txt')
do i = 1, niter1
WRITE(1,102) hist1(i)
end do
else
write(*,*) "Error: No sign change or failure on [-4,-2]."
end if
CLOSE(1)
CALL flush(1)
call bisect_hist(hist2, a2, b2, tol, niter2)
if (niter2 >= 1) then
write(*,'(A,1X,F12.7,1X,A,1X,ES12.4)') "Root on [2,4]:", hist1(niter2-1), "f(x*)=", f(hist2(niter2-1))
write(*,'(A,1X,I0)') "Iterations:", niter2
else
write(*,*) "Error: No sign change or failure on [2,4]."
end if
102 FORMAT(F0.17) ! Floating point format specifier, with 17 decimal places; the 0 here means that processor selects the smallest positive field width necessary
end program main
! ---------- external subroutine ----------
subroutine bisect_hist(hist, a, b, tol, it)
implicit none
real, dimension(:), intent(inout) :: hist
real, intent(inout) :: a, b
real, intent(in) :: tol
integer, intent(out) :: it
real :: mid
it = 0
write(*,'(A," [a,b] = [",F0.6,", ",F0.6,"]")') "Execution for interval:", a, b
if ( f(a) * f(b) >= 0.0 ) then
write(*,'(A," [a,b] = [",F0.6,", ",F0.6,"]")') &
"f(a)f(b)<0 not satisfied! f has no root in", a, b
it = -1
return
end if
do while ( abs(b - a) > tol )
mid = 0.5 * (a + b)
it = it + 1
hist(it) = mid
if ( f(a) * f(mid) < 0.0 ) then
b = mid
else
a = mid
end if
end do
end subroutine bisect_hist
! ---------- external function ----------
real function f(x)
implicit none
real, intent(in) :: x
f = exp(-x) * ( 3.2 * sin(x) - 0.5 * cos(x) ) - 3.0
end function f
And the separate params.f90 module:
module params
implicit none
real :: a1 = -4.0, b1 = -2.0
real :: a2 = 2.0, b2 = 4.0
integer, parameter :: n = 50
end module params
These can be compiled together with:
gfortran -ffree-form -c params.f90 bisection.f90
And then linked with:
gfortran params.o bisection.o -o bisection
which produces the executable bisection that can be executed with:
./bisection
Common mistakes#
Here is a list of common mistakes that a few people made:
Not using the absolute value to compute the length of the interval \([a,b]\). What if the values were complex? Distances (or lengths) are non-negative numbers by definition.
Using
while ((b - a) / 2.0 > tol)as a stopping criterion, effectively doubling the toleranceUsing the midpoint of the interval
abs(f(mid_point)) > tolto check for the stopping criterionNot checking the necessary condition for the function to have a root in the interval. That is, not checking that it has opposite signs at the endpoints
f(a)*f(b)<0
Part 4 and EC#
For Part 4 and the Extra Credit question, here are the Julia codes:
using Plots
default(linewidth=4, legendfontsize=12)
#### Part 4: Conditioning
f(x) = exp(-x) * ( 3.2 * sin(x) - 0.5 * cos(x) ) - 3.0
plot(f, label = "f(x)")
fprime(x) = -exp(-x)*(3.2*sin(x)- 0.5*cos(x)) + exp(-x)*(3.2*cos(x)+0.5*sin(x))
kappa(x) = abs(fprime(x)) * abs(x) / abs(f(x))
plot(x -> kappa(x), ylims=(0,20), label = "kappa for f(x)")
# savefig("kappa.png") # need to comment for published website
We can see how the function is ill-conditioned at \(x=-3\), where the root \(f(x^\star)=0\) lies, and well-conditioned in the \([2,4]\) interval.
### EC: Convergence rate
hist = zeros(15)
# open output.txt data and store values in the array called hist
open("../fortran_programs/module5-7_hw2_review/output.txt","r") do f
line = 1
while ! eof(f)
l = readline(f)
hist[line] = parse.(Float64,l)
line += 1
end
end
r = hist[end]
scatter( abs.(hist .- r) .+ eps(Float64), yscale=:log10)
ks = 1:length(hist)
plot!(ks, 7 * (.5 .^ ks))
# savefig("error_plot.png") # need to comment for published website
From the graph, we can see that the bisection algorithm is linearly convergent. In fact, we saw in class that it is \(r\)-linearly convergent, i.e., there are some errors in the sequence that are not strictly decreasing.
2. Submission expectations#
In general, when you receive an Assignment via a GitHub Classroom link, you want to clone your assignment repository, by doing
git clone your_assignment_repository_url
You can also work on a back-up repository or directory of your choice if you want to, for your scrap work, but you have to clone the assignment repository and submit your work there to be considered for submission and grading.
As soon as you clone your Assignment repository, move to that repository
cd your_assignment_repository
Create a new feature branch and switch to that. You can do this in two ways:
git checkout -b name_of_your_branchgit branch name_of_your_branchand thengit checkout name_of_your_branch
Do NOT work directly off
mainYou can work on your feature branch as much as you like and create repeated incremental snapshots of your work via
git commit. Always remember to use meaningful commit messages to remind yourself (and others) about your work in that moment in time. In a terminal you can simply do this by
git commit -m "Your commit message"
You can also write multi-line more detailed commit messages if you want. Just simply separate them with a space, and repeat the -m option, as in:
git commit -m "Your commit message" -m "Your more detailed message on a new line"
When you are satisfied with your committed work, you can push it to your working branch via:
git push origin your_branch_name
If it is the first time you are doing this, git will automatically tell you that you can open a Pull Request with your changes. Just CTRL-click on the URL that git shows you in the terminal and you will be sent to your Pull Request web interface.
Any successive changes that you want to push to your branch, they will be automatically reflected on the open PR.
Only changes made within the deadline (including the lateness window) will be graded.
Remember not to attempt to close or merge your PR without any Reviewer (in this case your instructor) approval.
Always remember to double check the
File changedtab in your PR. If you see files that should not belong there (e.g., files automatically created by your IDE or virtual environment files) remove them.If you are using an IDE that automatically creates hidden project files that you might inadvertently push to your branch, it is always a good practice to use a
.gitignorefile that specify which files you do not want to be tracked bygit, and therefore, pushed to your branch. Recall that we covered this in our first lecture.
2.1 Assignments in Julia#
If you are asked to submit your code in Julia, as a
.jlfile (not a Jupyter Notebook.ipynbfile), recall what we have covered during our first lecture.You are required to submit not only your Julia code (i.e., the file with the
.jlextension), but also theProject.tomlandManifest.tomlfiles that are created when you instantiate your environment. These files will be automatically populated if you need to use any Julia package that is available in the Julia registry, meaning, any package that you add by doing
]add PackageName
or in the Julia REPL:
using Pkg Pkg.add("PackageName")
These files are fundamental for your reviewer to be able to reproduce exactly your environment, run your code, and reproduce your results.