Untitled
Exercise 5 (Getting started)
a) Compute the expressions 127/13+pi, exp(3)*13i and exp(pi*i).
) Let C =
0.5
2
1
und B =
2 3 2
−0.25 0 −1
3/5 −3/5 0
. Determine C ∗ CT and CT ∗ B ∗ C.
c) Determine 23+10^ XXXXXXXXXX^37 by hand and numerically. Explain the difference.
Exercise 6 (Using commands)
Let C, B be as in the previous exercise.
a) Compute the norm and the mean value of C.
) Compute the determinant and the inverse of B.
c) Determine the eigenvalues and the eigenvectors of B.
Hint: For unfamiliar commands use help or check on the web.
Exercise 7 (Extracting matrix entries)
Generate a 5 ×5 matrix A using the command magic and determine the following expres-
sions:
A(:,4) A(2,:) A(1:2:5,:)
A(2:4,2:4) A(3,5) A([1,5],[1,5])
Exercise 8 (Operators)
Let the variables a = 5, b = 2 and c = −1 be given. Firstly, determine the results of the
following expressions by hand; secondly, check them using Matlab.
a < 10 & a > 0 a < 4 | ( c <= 2 &
2 ) c ~= -1
== a + 3*c abs(c) -
2 > 0
Exercise 9 (Function functions)
Let the following piecewise continuous function be given:
f(x) =
− sin(x) x < 0
x2 0 ≤ x ≤ 1
1
x
x > 1
a) Write a function file that computes f(x) (signature: function y = f stk(x)). Ensure
that a vector valued call of the function is possible.
) Plot the function on the interval [−π
2
, π] using plot.
1
Matlab Computer Exercises (II)
Exercise 10 (Threedimensional plots)
a) Plot the function x = cos(t), y = sin(t), z = cos(2t) on the interval [0, 2π] using the
command plot3.
) Create a surface by means of [X,Y,Z] = peaks(30) and plot it with mesh, plot3 and
surf.
c) Visualize the function
g(x, y) = log
(
1
√
x2 + y2
)
on the domain [−1
2
, 1
2
] × [−1
2
, 1
2
] using mesh.
Hint: The co
esponding matrices can be constructed employing linspace and
meshgrid.
Exercise 11 (Programming exercise)
Write a program that randomly draws m = 10 regular polygons with n = 3, . . . , 30 nodes.
The nodes of each single polygon must lie on a circle with radius 1. The midpoints of the
polygons shall be distributed randomly in the domain [0, 10]× [0, 10], and shall be drawn
in a randomly chosen color.
Hint: Find out about the function patch by using the Matlab help.
2
Exercise 12 (Explicit one step methods)
a) Implement the improved Euler method in a m-file with the signature
function [t,y] = IMPR EULER(FUNC,t0,T,y0,n,var).
Make sure that your code can solve systems of ODEs where the unknown y ∈ Rd is a
column vector.
) Implement the linearized pendulum ODE which is given by
y′
1
(t) = y2(t),
y′
2
(t) = −
g
l
y1(t).
The file should have the signature [f] = LINEARPENDULUMODE(t,y,var) where va
contains the positive constants g and l.
c) Copy the file EULER.m from the last project into your cu
ent directory. Download the
script pendulum comp.m from the web page and run it. Explain the results.
Exercise 13 (Fixed point iteration)
We want to determine the root of a nonlinear function g(x) by means of a fixed point
iteration.
a) Let g(x) = x3 + 4x2 − 10. Show that any solution x̄ with g(x̄) = 0 is a fixed point of
each of the following functions:
3
(i) φ1(x) = x − x
3 − 4x2 + 10
(ii) φ2(x) =
1
2
(10 − x3)
1
2
(iii) φ3(x) = x −
x3+4x2−10
3x2+8x
Hint: A point x̄ is a fixed point of a function φ if φ(x̄) = x̄.
) Implement the fixed point iteration method in a m-file with the signature
function [fp,iter] = fixed point(PHI,x0,TOL,Nmax).
The first input argument is the function handle PHI for the fixed point iteration;
further, x0 is the initial approximation to the location of the fixed point, TOL is the
absolute e
or tolerance, and Nmax the maximum number of iterations to be performed.
The function should return the approximate value of the fixed point fp and the numbe
of iterations iter.
c) Test your implementation for the function g from a), which has a solution in [1,2].
Choose TOL = 10−5, x0 = 1.5 and Nmax = 15. Run the iteration with φ1, φ2 and φ3
and compare the results.
4
Exercise 14 (Classical Runge–Kutta method)
a) Implement the classical Runge–Kutta method with the Butcher table
0
1
2
1
2
1
2
0 1
2
1 0 0 1
1
6
1
3
1
3
1
6
in a m-file with the signature
function [t,y] = RK4class(FUNC,t0,T,y0,n,var).
As in Exercise 12, the function should return the vector of time steps t = (tj )j=0,...,n and
the co
esponding approximate solutions y = (yj)j=0,...,n which can be vector-valued.
The input parameters are the right hand side function FUNC, the time interval [t0, T ],
the starting solution y0, the number of steps n and the additional parameter var that
can pass further variables to the ODE.
) Modify the script pendulum comp.m from Exercise 12 such that the RK4 method of a)
is included in the comparison. Compare the numerical convergence orders.
c) Change the script of b) such that the scalar IVP y′(t) = t0.5, y(0) = 0, is solved.
Compute the e
or at time T = 1 where the exact solution is given by y(1) = 3
2
. Do
the methods still yield the same numerical convergence order as in b)?
Exercise 15 (Stability of RK methods)
a) In the lecture, the stability function g(z), z = hλ, of a general Runge–Kutta method
is given by the formula
g(z) = 1 + zγT (I − zB)−11.
Determine the stability functions for each of the following Runge–Kutta methods:
(i) Explicit Euler (ii) Improved Eule
(iii) Implicit Euler (iv) Trapezoidal Rule
Plot the absolute value |g(z)| for the above methods on the interval [−10, 0].
) Determine the minimal value zmin < 0 such that the condition |g(z)| < 1 is satisfied fo
all z ∈ (zmin, 0). For λ = −100, what is the largest step size hmax that can be chosen
such that the method still gives quantitatively co
ect results?
5
Exercise 16 (Step size control)
In this exercise, the classical RK4 method of Exercise 14 is extended to a RK4(3) method
with step size control.
a) Write a function [ynew,est,k1new] = RK43(FUNC,told,yold,h,k1,var) that per-
forms one step of the FSAL (First Same As Last) Runge–Kutta scheme with the
Butcher table
0
1
2
1
2
1
2
0 1
2
1 0 0 1
1 1
6
1
3
1
3
1
6
p = 4 1
6
1
3
1
3
1
6
0
p̂ = 3 1
6
1
3
1
3
0 1
6
The input variable FUNC contains the function handle for the right hand side of the
ODE. The other input arguments are told= tj , the last approximation yold= y
j, the
cu
ent step size h= hj and the value k1= k
j
1
computed in the last step. var contains
the parameters for the right hand side FUNC.
Inside the function, the fourth order approximation yj+1 as well as the third orde
approximation ŷj+1 are computed according to the formulas
yj+1 = yj + hj
5
∑
l=1
γlkl, ŷ
j+1 = yj + hj
5
∑
l=1
γ̂lkl,
with γ̂l given in the last line of the above Butcher table.
The output arguments are the new approximate solution ynew= yj+1, the e
or indi-
cator est= ‖yj+1 − ŷj+1‖ as well as the FSAL transfer k1new= kj
5
= kj+1
1
.
) Implement a function
[t,y] = SOLVE RK43(FUNC,t0,T,y0,var,tol,h0,param)
which solves the ODE with the right hand side FUNC by means of the routine of b). The
input arguments are the time interval [t0, T], the initial value y0, the parameters
var for the evaluation of FUNC, the relative tolerance tol, the initial step size h0 as well
as the specifications param= (hmin, hmax, αmin, αmax, β). These values are necessary fo
the adaptive calculation of the new step size hnew according to the following criteria:
qold =
est
tol · hold
,
α = min
(
αmax, max
(
αmin, q
−1/p
old
))
,
hnew = min
(
hmax, T − t, max(hmin, α · β · hold)
)
,
6
where p is the order of the method.
If qold ≤ 1 or hold = hmin holds, the approximate solution y
j+1 computed by RK43 is
accepted and the actual time t is increased to t + hold; otherwise, the value of y
j+1
is discarded and the calculation of RK43 is repeated with the new (smaller) step size
hnew. In either case, the new step size hnew is calculated and used for the next step.
The function should return the vector of time steps t and the co
esponding approxi-
mate solutions y which can be vector-valued.
c) Download the file threebody.zip from the web page, unzip it and start it by the
command threebody. Test your implementation by computing the o
it “Tuft” with
a tolerance of 10−4.
7