Scientific Computing 2018/19

Office hours:

Homeworks - lab problems

References

1. P. Krzyżanowski, Obliczenia inżynierskie i naukowe. Skuteczne, szybkie, efektowne, Wydawnictwo Naukowe PWN, 2011
2. P. Krzyżanowski, Obliczenia naukowe. Html lecture notes (in Polish). Pdf version is also available there. University of Warsaw. 2010.

• Lab - linux basics: cd, ls, rm, cp, pwd, mv, tar, mkdir, rmdir, chmod, head, tail, more, cat, man. Basics of octave - octave as a scientific calculator, semicolon : how to create a matrix, saving/loading from/to files in different formats; basic matrix/vector operations - opt. solving linear equations and LLSP etc Pr 1 Create matrices 3x4 A and 3x5 B - then the matrix 3x8 C with 3 first columns are A and next B. Create D 3x3 as a main minor of C oi.e. from C(1,1) to C(3,3). Change the order of columns of D. Compute sin on elements of D. Write D to a file (binary or ascii) - change D(1,1) to -100 and load this new matrix to octave as DD. Compute 1st norm of DD. Pr 2: Compute a discrete max norm of (sin(x))^2 on [0,1] (on 10000 points) Create a text file with any 3x3 matrix in a proper octave format and load it to octave. Plot a graph of sin(x) on [0,4]. Do the same for x*(sin(x)+2). Compute an approximation of derivative by central difference on [0,4] with h=4/N N=100 i.e. compute Dcf(xk)=0.5(f(xk+h)-f(xk-h))/h for xk=4*k*h k=1,..,N (equidistant points on [0,4]. Can it be done without a loop? Compare with the values of real derivative at those points i.e. er(xk;h)=|Dcf(xk)-f'(xk)|. Check the order of approximation i.e. take once N=100 then N=200 and plot er(xk;h)/er(xk;h/2).
  Lab - scripts: example matbasic.m write on command line of octave'a: matbasic.
• Lab - functions, m-files, basic linear algebra, solving the systems, LLSP, norms, matrix factorizations,: Cholesky, LU, QR, eigenvalue problems etc
  • Write a m-file with a function which has two parameters x and a and returns f(x;a)=x^a and its derivative as a 2nd returned value; default value of a=1.
  • Solve any linear 2x2 system e.g. A=[1,2;3,4] f=[1;3], x=A\f - test if residual vector f-A*x is zero, and if we get the right solution (if we know it) here the solution =[1;0]
  • Solve any simple regular LLSP - e.g. A=[1,1;1,2;3,4] f- the 1st column of A - x=A\f - is the computed solution equal to the 1st unit vector ie. to [1;0]?
  • A simple problem - polynomial interpolation; adn linear regression as a LLSP
    • using vander() and polyval() and backslash - solve linear regression problem for givem x_k,y_k find a,b such that \sum_k |ax_k+b -y_k|^2 is minimal. Compare with polyfit(). Write a function which for two parametes being the vecors x and y (of the same length) returns the vector wit hthe parameters [a;b] of the linear regression problem.
      More demanding version : The function should check if the problem has a unique solution (up to fl calculations naturally)? How to do it withou using rank(), svd(), or null() functions? HINT: perhaps we can use qr() instead of backslash?
    • find Lagrange interpolating polynomial on equidistant points on [-5,5] for f(x)=1(1+x*x) - compute discrete max norm of the error on [-5,5] and at interpolation nodes Test for 5,10,20 and 30 nodes.
    • do same as in the previous problem but for Chebyshev points (linearly scaled roots of T_n(x)=cos (n*arc cos(x))).
    • Not always direct solvers gives a good solution. Take N=2+5*k equidistant nodes on [0,1] and repeat computing Lagrange interpolator p for any function e.g. f(x)=sin(3x) for k=0,1,2,3,...; check if f=p at interpolation nodes i..e compute max_k|f(x(k)-p(x(k))| where (x(k))_k interpolation nodes.
  • Take a symmetric matrix 2x2 and compute its eigenvalues and eigenvectors. Check if AX=XD X- a matrix having the eigenvectors as columns - D a diagonal matrix with with respective eigenvalues on its diagonal. Test for A=[2,1;1,2], then D=diag([3;1]) and X=[1,1;1,-1]. Compare with the results of eig(). Repeat tests for different random symmetric matrices.
  • SVD: compute svd() for a random matrices Check if the obtained values are OK.
  • Solve a LLSP problem by svd() for a random matrix 5x4 and 4x5. Compare with backslash results.
  • Compute using svd() condition of A for a random A 5x5; compute the second matrix norm of A for a random A 5x4 etc Compute the orhonormal basis of ker(A) and R(A). Compare with the results obtained by null() and orth().
  
  Some solutions
  lab2.m - a few solutions of LLSP and linear systems (interpolation), linear regression
  

  Home problem no 1.(graded) Due Tuesday March 19, 2019

  Using octave functions (no loops!) write a m-file ltrigp.m with the function [sol,info]=ltrigp(x,y) solving the following problem: for a given two vectors x and y (of the same length) compute [the vector [a;b;c] such that a,b such that \sum_k |a+ b*cos(x_k)+ c*sin(x_k)-y_k|^2 is minimal. In case when a solution exists but is not unique the function should return the solution of the minimal 2nd norm. Thus the function should have 2 parameters: the vectors x,y and should return sol - the vector with the solution and info the code: 0 - there is a unique solution; 1 - there is no unique solution; 2 - the entry vectors have different length. In case of info<>0 the function should print a message to the display.
  One has to show me the working code in the lab on Tuesday, 1215pm-145pm, March 12 or 19, 2019.
  HINT: you can rewrote the linear regression function which we should write in the lab (I hope...) It is also enough to check numerically if the revelant matrix is singular or not-singular (or rather of max rank or not max rank).
• Lab - chapter 5 breeding of animals i.e. applied linear algebra. Link to pdf (unfortunately only in Polish).
  Prob (NOT graded).
  Given 3 matrices X,Z,A=A'>0 and parameter k>=0 write 4 solvers for Gs=r with G=[X'*X,X'*Z,Z'*X,k*A^(-1)] and r=[X'*y;Z'*y]
  1. a direct one - create matrix (using inv(A)) and the rhs vector - use backslash
  2. solving an equivalent LLSP using backslashem - with a block upper triangular matrix [X,Z;0,S] and the RHS vector [y;0] for S=S'=sqrt(k)*(sqrt(A)^(-1) - the square root of a matrix can be computed by sqrtm() or by solving an eigenproblem
  3. solving an equivalent LLSP using backslashem - with a block upper triangular matrix [X,Z;0,S] and the RHS vector [y;0] for S=sqrt(k)*L^(-1) -with L the Choleskiego factor of A (R=chol(A);L=R')
  4. solving an equivalent LLSP using backslashem getting a vector s=[b;g] - with a block upper triangular matrix a [X,Z*L;0,S] and the RHS vector [y;0] for S=sqrt(k)*Id -with L the Choleskiego factor of A (R=chol(A);L=R')- then 2nd term of the solution is a=L*g i.e. g=s(3:end)
  Compare the time and if we get the same solutions for example 5.1, m=5,n=8,k=2 matrices are in the lecture notes -one can also use: load rozd5data.txt - the file with the data of this example 5.1.
  Some solutions
  lab3.m - solutions of the breeding problem
• Lab: Numerical linear algebra cont. Simple plots, boolean operators. Sparse matrices in octave. boolean operators (&& and, || - or, !-not, != - different, == - equality itp), in particular boolean operator element by element
  Solve nonlinear equations
  1. Prob create a 3diagonal matrix L with 2 on the main diagonal and -1 on the sub- and super diagonals using different methods: loop,loop + blocks, function diag(), sparse() etc. Compare time: - tic and toc.
  2. Prob Compute an inverse of L. Solve the system Lx=f in 3 different ways: 1/ using backslas, 2/ using backslas wit hL as a sparse, 3/ multiplying f by the inverse of L - compare time (inverse of L is precomputed so measure only the time of multiplying by it)
  3. Prob: implement the function hat without loops using only vector operations (f(x)=f(-x);f(x)=0 outside [-1,1], f(x) =1-x on [0,1]).
  Some solutions
  lab4.m - some solutions - graph of max(f,g) without loops, creating sparse matrices - comparing times of soliving linear systems etc
• Lab: Numerical classics, solving nonlinear equations etc global and static (persistent a=0) variables, (fsolve and fzeros) for simple equations x^2-2=0 etc exp(x)+x=0 . Global variables can be useful....
  • Solve x*(1+0.5sin(x))=0; (x-1)^2=0;
  • using fsolve plot the graph of implicit function: e.g. y(x) for y^2+x^2=1 y(0)=1 or y(0)=-1; another ex. compute y(x) for f(x,y)=x^3+y^3-y=0 near y(0)=1.
  • plot the graph of the inverse of 2x-cos(x) on [-2,2] or tan(x*x) on [0,60] or x^ on [0,2] etc
  
  Some solutions
  lab5.m - a script with solutions of some problems
• Lab Nonlinear problems cont. Conditionals, loops etc if else, switch, loops: while, do ... until, logicals (&& and, || - or, !-not, != - not equal == equality), in particular x||y for x,y, vectors/matrices
  • Repeat in 2d : on the mesh 10x10 on [0,1]x[0,1] compute the values of the inverse function of f([x;y])= [10x+y*y; 10y+x*y]
  • A bit more serious project we want to plot the inverse of a given difficult function - see formulas in Chapter 3 in lect. notes.

  • Home problem no 2.(graded) Due Tuesday April 16, 2019

    Using octave functions write a m-file invgraph.m with the function [if,y,info]=ltrigp(fun,a,b,N) computing the approximated values of the inverse of a given univariate function: f:[a,b] -- >R; over N equidistant points over the interval [f(a),f(b)] (or [f(b),f(a)]).
    INPUT:
    1. fun- a function handle to an univariate function f of the form y=f(x);
    2. a,b, the ends of interval
    3. N positive integer value;
    OUTPUT
    1. if the vector of approximated values of f^(-1) over the points in vector y
    2. y- the equadistant points over [f(a),f(b)] or [f(a),f(b)]
    3. info - exit code - 0 if the sucess i.e. the values can be computed; nonzero - if something wrong (e.g .f(a)=f(b) or the solver failed etc)
    Please test on some simple invertible functions like x^2 on[1,2], 1-x^2 on [1,2], sin(x) on [-1.5,1.5], 2x+sin(x) on [-1,2] - plot the graphs of f and inverse of f.
    One has to show me the working code in the lab on Tuesday, 1215pm-145pm, April 16, 2019 or earlier.
    HINT: you can use the solution of the problem from Chapter 3 of lecture notes.
  • Prob : implement without loops a hat function (f(x)=f(-x);f(x)=0 outside[-1,1], f(x) =1-x on [0,1]) or 1_[0,1/2](x) -the characteristic function of [0,1/2].
  • Quadratures - numerical integration ( quad(), ), integral with singularities
  • Prob compute integrals of a few simple functions e.g. sin(x), polynomials: x^2, x^p etc integrals over unbounded intervals e.g. \int_{-\infty}^{+\infty} exp(-x^2)dx or integrals of a function with a singularity e.g. \int_{-1}^1 \sqrt{|x|}dx (no differentiable at zero) or \int_{-1}^1 |x|^pdx dla -1 <= p <= 1.
    
  • Prob Plot the graph over [-20,20] of F(x)=\int_{-\infty}^x exp(-t*t/2) dt. Check if lim_{x tends to \infty}F(x) == \sqrt(2\pi) and lim _{x-- > -\infty}F(x)=0. Solve F(x)=1/2. How to write an octave function computing the integral L^2 norm over [a,b] - 3 parameters: a function handle, and a and b
  • Prob Write a similar function computing L^p (p-2nd parameter).
  • Prob : Write a function computing \int_0^1f(x)x^k dx for y=f(x) and then using also hilb function solve a polynomial approximation problem in for a hut function or 1_[0,1/2]. (i.e. solve the system with the Gramm matrix in L^2(0,1) for (1,x,x^2,..,x^m) and the rhs: f_k=\int_0^1 f(x)x^k dx k=0,..,m.
  • Prob : Compute a 2d integral \int_0^1\int_0^1 sin(x*y) dx dy or \int_([0,1]^2) sin(x+y^2)*(x-3y) dx dy or the 2d volume of the elipsis 2x^2+3y^2-1=0 (y(x) is given implicitely... ) - you can use either dblquad() or just quad() twice...
  • Prob : Compute a 3d integral \int_0^1\int_0^1\int_0^1 sin(x*y*z) dx dy or \int_([0,1]^2) sin(x+y^2-z)*(x-3y*z) dx dy or the 3d volume of elipsoid 9x^2+4y^2+z^2-1=0 (z(x,y) is given implicitely... ) - you can use either triplequad() or just quad() three times... compare with the analytical value (it can be easily computed by a substitution F(x,y,z)=(3x,2y,z) which maps it to the unit ball)
  
  lab6.m - the solutions of problem from chapter 3 in lecture notes and some simple examples of compting integrals- simple and double...
• Lab. Ordinary Differential Equations- Initial Value Problems
  Solving initial value problems for ODEe : main solver : lsode. - scalar and multivariate cases
  • Prob 1 Solve a few simple ODEs classics e.g. y'=-a*y;y(0)=-y on [0,1] or[-10,0]. Blow up equation: y'=y*y; y(0)=1 on [0,1/2] or [0,3]. Linear penduulum equation: y''=-y with y(0)=1;y'(0)=0 non[0,2*pi](solution: cos(x)) or a nonlinear one: y''=-sin(y) with the same initial data (is the solution periodic?) . Solve (lsode) stiff system using bdf and adams - copare time (lsode_options) e.g. : x'=998*x+1998*y;y'=-999*x-1999*y ; x(0)=y(0)=1 on [0,10], [0,100] etc Compute eigenvalues of the matrix [998 1998;-999, -1999].
    or the stiff system of chemical reaction of Richardson: x'=-0.04x+10^4*y*z;y'=0.04*x-10^4*y*z-3*10^7y*y; z'=3*10^7*y*y on [0,40] or [0,10^k] k=5,6 etc -try to solve it using standard lsode options and options for explicit schemes. Compare time.
  • Prob Plot the graph of F(s)=y(1,s) with F(s)=y(1,s) where y(t,s) solution of: y'=f(y,t), y(0)=s for s in [0,1], [-1,1] for f=y,10y,-10y 0.1*y*y itp
  • Prob Consider IVP: y'=y*(1-y)*sin(y*y+x); y(0)=s. Let y(x,s) be its solution for s\in[0,1].
    • Compute the integral \int_0^1 y(1,s) ds.
    • Find z such that \int_0^z y(1,s) ds=0.2 and find s0 such that y(1,s0)=0.5.
    • Plot the graph of s|-->y(1,s).
    • Find s0: y(1,s0)=0.5 by solving a revelant nonlinear equation. I.e more generally speaking how to compute s for a given z\in[0,1] such that y(1,s)=z (if we define function F(s)=y(1,s) then z is the value of the inverse of F at z). Is it possible to find z with one call to lsode() i.e. without solving any equations?
  • Prob Blow up. Solve two problems u'=100*u u(0)=1 on [0,10] and [0,100] and u'=u*u u(0)=1 on [0,10].
  • Prob 200,00 Solve boundary value problem: : u''+u*u=f(t) u(0)=u(1)=0 for different RHS function f e.g. f=1, f=exp(t) etc Plot the graph of the solution. Find min and max of the solution. By the way learn how to plot graphs in semilogs scales (semilogy, semilogx), parameters of plot etc How to do it?
  
  lab7.m - some simple solutions of solving IVP with lsode().
• Lab. An Earthquake model of a building with stiff floors - with elastic connections between them and the zero floor with the ground, assume that the i+1 th floor of the mass m_{i+1} which moves relatively to the i-th floor: x_{i+1}-x_{i} then a restoration force k_{i+1}(x_{i+1}-x_i) (proportional to the relative shift) acts on the both floors i.e. on the i-th acts 2 forces:

  m_i*x_i''=k_{i+1}(x_{i+1}-x_i) -k_i(x_{i}-x_{i-1}), forces on the i-th floor, i>0

  on the zero floor (ground floor) acts also analogous restorative force proportional to the shift relative to the ground: i.e.

  -k_0x_0

  and the restorative force between zero floor and 1st floor: i.e. m_0x_0''=k_1(x_1-x_0) -k_0x_0. Naturally, on the last floor n-th acts only one restorative force... m_n*x_n''=-k_n(x_n-x_{n-1}).
  During than earthquake additionaly the force

  f(t)=G*sin(gamma t)

  acts on the ground floor (level 0) e.g .gamma=3.
  we have ODE: MX'=SX+F for M diagonal matrix n+1 x n+1 (masses of floors on the diagonal), matrix S a 3-diagonal symmetric matrix such that on the main diagonal we have -k_{i+1}-k_i (k_{n+1}=0) sub- super- diagonal = (k_1,..,k_n)
  F=(f(t),0...,0)' z f_1(t)=Gsin(gamma*t) - G<>0 during earthquake (in practise f should of the form f(t)=G*sin(gamma*t) exp(-a|t-t_0|) a - a positive large constant. Problems:
  • Write ODE: M*X''=S*X+ C*X' + F(t). (F =(f(t),0..,0)^T, C - diagonal negative matrix - air/friction resistance. Write functions computing matrices: A=M^{-1}S, M^{-1}C and the vector function M^{-1}F
  • For 2 storey building m_k=m=10000kg z k_l=k=5000kg/s^2 compute eigenvalues of M^{-1}S : \lambda_l - and next the frequancies of the building defined as \omega_l=\sqrt{-\lambda_l} and periods T_l=2\pi/\omega_l. Are the periods in the range [2,3]? (a typical range of an earthquake)
  • Simulate 2 floor bldg with some initial nonzero (small) values without earthquake (so as the building was moving after an earthquake). Plot the graphs of x_i.
  • Do the same but for for 10 floor building( all floors has the mass = 10000kg) and k_l=k=5000kg. Are the periods in [2,3]? Simulate the solutions for t in [0,200]. Repeat adding damping: C=-alfa*Id for alfa=0.1,1,10
  • Compute frequencies and periods of the pyramid building i.e. such that m_1= 10000kg and m_{i+1}=0.8*m_i (k_i = 5000kg/s^2) Compare with a standard 10 floor building with m_i=m=10000kg and k_i=k=5000gh/s^2 - Compute frequencies and periods - are periods in the range of 2-3 seconds (typical for an earthquake)? . Simulate an earthquake for gamma=3 and G=10,20,40,80cm. Plot the graph of the highest floor - compute max_i |x_i(t)| Add damping alfa=0.1,1,10. Change G = G(t) e.g. G_{max}exp(-t^2) or G constant for 10secs and then rapidly decreasing
  • Compute the same for a upside down pyramid i.e. m_i=0.8*m_{i+1}....

  • Home problem no 3.(graded) Due Tuesday May 21, 2019

    Consider a 10th story building with a swimming pool in the ground floor. We want to solve numerically a model of an earthquake acting on this building (model above) without damping $C=0$ The mass of a store no. 1-9 m(j)=10000kg for j=1,...,9. and of the first one: m(0)=M=30000kg, constants k(0)=10000kg/s^2 and k(i)=5000kg/s^2 i=1,..,9 - the earthquake force acting on the ground floor: 5000*cos(3*t) kg/s^2 times cm. Goal: compute the frequencies and periods of this building. Check if periods are in [2,3]. Plot the graphs (in one graphical window (figure)) of the displacements of 1st and the last floors (x1(t) and x10(t)) on $[0,T] for T=15. We assume that x(0)=x'(0)=0. Compute the amplitudes of the shifts of all floors (for each independently) i.e. compute (max_t |x(k)(t)|)_k k=1,..,10 and plot the graph of F(t)= max_k|x(k)(t)|.
    
• C language - simple codes (one can use C++ or fortran)
  • Prob Write a simple C program which does something e,g, let it print a message "Hello world" or asks for two number - and when user provides them the program prints a sum on the screen. Compile it and execute it.
  • Write a code in C/C++ which sums the harmonic series:\sum 1/n^p; p>1 using different loops: for(){ }, while(){ }, do{ } while(), option - p should be given interactively from the keyboard by an user
  • Binary and ascii files in C. Read help to (internet or in unix: man fopen )fopen(), fclose(), fwrite(), fread(), fprintf(), fscanf(), fgetchar() feof()etc Write a code that will write a vector (double array)to a file - binary or ascii (text) file then another code that will read it from the file. Test for: [1 2 3 4 -9.45].
  • Write a program that will read to memory a matrix saved from octave in ascii or text format (homework, not graded)
  • Write a macro which enables access to a matrix entries kept columnwise (or rowwise) in a single array (so as a vector). Then write a function multiplying a vector x by A or A^T (A kept columnwise in a 1d array)
  • Prob Memory allocation. allocate and then free dynamically vectors float,double etc. Write a a C function generating a vector sin(x) for a given vector x (e.g. equidistant points on [a,b]) - parameters : a pointer to x (or ends of interval - easier version)) and pointer to returned vector with sin(x) - C function should allocate memory and write to that memory sin(x), optionally export this vector to binary file (fwrite()) and then read it back to memory to another memory address.
  • Precompiler macros like functions - #define #if etc Modify the code from the previous function in such a way that depending on the precompiler options it will work either for doubles or floats (single precision).
  
  lab9-C.zip - zip file with solutions of a few problems (lab?a.c reading parameters when the executable program is run); unpack it with unzip lab9-C.zip command
  Homework (not graded.. just for yourself)) (1) write a C code that reads matrices from a text file written in ascii octave format to a double indexed float array in C. Then show the matrix on the computer display, i.e. save a matrix in text format in octave and then display by your C programm - check if it is OK. (2)Write another C code which writes a matrix from C to text file in a default text format of octave. Read it using octave - check if it works. (3)Try another wrtitng to other octave formats (ascii files) - or to binary octave format - (difficult - one has to find details of octave binary format)
• Using fortran libs i.e. e.g. BLAS or LAPACK in C/C++
  • Prob write a macro IJ(i,j) to a matrix MXN which is held in memory columnwise as a single vector of length M*N, and a functions that using that macro can show this matrix on the display, (if you are fast you can also write functions such that : adds two matrices (in that format), multiplies two matrices, scales a matrix, creates matrix of given dimension (allocates memory), destroys matrix (frees memory) etc )
  • Prob create a new type which is a structure with two fields: 1st integer (dimension of a vector), and 2nd real pointer (data)- which represents a vector - the na functions creating vector (allocating memory to data pointer), destroying, scaling a vector, adding two vectors etc. (homewrok - not graded)
  • Prob use a fortran BLAS function e.g. dnrm2.f- computes a 2nd norm of a vector - or - dscal.f - scales a vector by a scalar - in C code. Then write your header file : (in C code: #include "mojblas.h") with the proper headers of this functions (in order to use them in your C code)
    In general we have to compile fortran libs/functions using lib: gfortran.g and BLAS using blas.a, e.g. gcc ourccode.c -lgfortran -lblas or if we want to compile using fortran e.g. function dscal.f then we can write: gcc ourccode.c dscal.f -lgfortran (-lm - in case of dnrm2.f - this funcion needs sqrt() from math lib) - now on some modern systems we can omit -lm or -lgfortran; they are added automatically
  • Prob Using dnrm2.f (computeds the 2nd norm of a vecotr) and ddot.f - computes a dot produt of two vectors - write a function computing the cosine between two vectors i.e .cos(x,y)=x'*y/(||x|| ||y||) here ||x|| - the 2nd norm of x.
  • Prob Using a macro: IJ(i,j,M) returning the index of the element (i,j) in a matrix MXN kept in a single index array columnwise M in C and BLAS DGEMV(): dgemv.f compute Ax for A=[1 2;-1 1] and x= [0;1], [1;0], [1;1].
  • Prob : Write your own header file myblas.h with a prototype of BLAS DNRM2() and other BLAS functions and compile your program calling this header file using a precompiler command: #include "myblas.h". Naturally, remove any BLAS function prototype from the C source file.
  • Prob : Solve system Ax=f for A=[ 2 -2; 1 1] and f=[0; 2] (rozw [ 1;1]) using LAPACK function DGESV(): dgesv.f
    Compute the inverse of A=ones(6,6)+10*eye(6) and solve AX=B for B s.t. 1st column of a known X e.g. X[k]=1.0 i.e. take B=A*X for a known X, then, solve the system AY=B - check if Y==X.
  • Prob : Using BLAS functions dgemv() dnrm2() daxpy() and lapack dgesv() - create Hilbert matrix H(1/(i+j-1))_{i,j}^N (this must be done in a double loop), then using dgesv() solve the system Hx=b for b being the 1st column of H (b(i)=1/i) - then using dgemv() compute r=H*x-b and using daxpy() e=x-s with s being the 1st unit vector - the solution (or just take: e[0]=x[0]-1) then using dnrm2() compute the 2nd norms of e and r and print them out on the display. Try N=5,10,15,20. (naturally you can keep r in b adn e in x, i.e. compute b=Hx-b and x=x-s).
  • Prob : Solve system Bx=f for B tridiagonal NxN (discrete Laplacian 1D) : 2 on the main diagonal and -1 using LAPACK function DPTSV(): dptsv.f , check if we get the right solution ( e1=[1;0;0;...;0]) using a right BLAS function.
  • Using LAPACK function dgtsv() solve the tridiagonal system system Ax=f with the tridiagonal matrix A with 7 on the main diagonal, -2 on superdiagonal, and -1 on subdiagonal for f=[7;-1;0;0;..;0] for N=15.Check if you get the right solution!
    dgtsv.f
  • Prob : Using a LAPACK function DGELS(): dgels.f for solving a regular LLSP Ax=f. Test on A=[1,1;1,2;1,3] and f=[2;3;4] (i.e. x=[1;1]) Solve the system MX2 Aa=y equivalent to linear regression i.e. x,y given vectors, we want to get a=arg\min_a \sum_k|a(1)+a(2)x(k)-y(k)|^2. Write a C/C++ function solving linear regression by DGELS(). Test for x=[1;...;10] and y=x+1 and different M=2,3,5,0,120. INPUT: integer M (dim of x,y), pointers to double x,y, OUTPUT: pointer to double a (solution of LLSP), pointer to integer info (code of results, zero if OK) (a=[1;1]).

  • Home problem no 4.(graded) Due on Tuesday June 4, 2019

    Write a function nellspsv() (in a separate file nellspsv.c) which solves a max rank LLSP with the matrix A and the rhs vector b (A m x n of maximal rank with m>=n) first forming the normal matrix B=A'*A and the vector g=A'*b using the respective BLAS functions (DGEMM() and DGEMV()) and then solving the resulting normal equation system by the lapack function DGESV(). (A' = A transpose) Then, write a code testing the function on some simple examples e.g. A=[1,1,1;1,2,3]' and b=[1,1,1]' (so the solution is x=[1;0]). During the grading I may ask for another LLSP (of the same size). Compile the both files (one with the test and nellspsv.c) and run it displaying the solution.
  
  testmacro.c - a simple program using a macro for accesing (i,j) element of a matrix stored in 1D array columnwise - there is a function printing the matrix to the display
  cosxy.c - a simple program with a function cosxy() computing the cos between two vectors x and y using two BLAS level 1 functions dnrm2() and ddot()
  tblas.c - a simple C code using BLAS dnrm2() function - and dynamically allocating memory (using malloc(); free() functions
• LAB: BVP solved by a shooting method and FDM method
  • Shooting method: solve BVP

    -u''+u=0; u(0)=u(T)=1

    Apply the shooting method, namely, solve IVP

    -u''+u=0; u(0)=1; u'(0)=s

    getting u(t,s) (practically numerically using lsode()). Find s0 s.t. u(t;s0)=1. It is enough to shoot twice i.e. take s2=1 and s1=0 and compute using lsode() u(t;0) and u(t;1). Then our s0=(1-u(T,0))/(u(T,1)-u(T,0)) (why?). Then compute u(T;s0) using lsode() and |u(T;s0)-1|. Why is the error growing with T?
  • Solve by the shooting method u''=(3/2)u^2 with u(0)=4 and u(1)=1. Define using lsode() the function F(s)=u(1;s)-1 with u(t;s) being the solution of IVP -u''=(3/2)u^2 with u(0)=4 and u'(0)=s. Solve the system F(s*)=0 using fzero() of fsolve() taking different starting s0 initial iterates for a nonlinear solver Plot the graph of F on [-100,0].
  • A simple eliptic problem:
    
    solve

    -u''+c*u =f
    u(a)=alpha, u(b)=beta

    (c - nonnegative constant) using 3point stensil FDM method. I.e. we are looking for uh(k) approximating u(x(k)) with x(k)=a+k*h h=(b-a)/N - uh(0)=alpha uh(N)=beta - solve the linear system AU=F with U=(uh(0),..,uh(N)) F=[alpha,f(x(1)),..,f(x(N-1)),beta] and A=A1+c*Id - A1 tridiagonal with A1(0,0)=A1(N,N)=1 A1(i,i)=2/h^2 and A1(i-1,i)=A1(i+1,i)=-1/h^2 i=1,..,N-1.
    • Create the matrix A with c=0 - display it for a small N e.g. N=6, do the same for F (taking e.g. alpha=beta=1, f=0, or f=sin(x), a=0,b=1) - solve the system using backslash - plot the graph of the FDM discrete solution
    • test the solver for a known solution i.e. take u=sin(x) then f=sin(x), alpha=sin(a) beta=sin(b) -test on a=0,b=1 taking N=20,40,80,60 - compute the pointwise error erN=|uh(k)- sin(x(k))| and its norm ||erN||_\infty and the discrete L^2 norm ||erN||_0^2= \sum_{k=0}^N h*|erN(k)|^2. Compute then the ratio erN/er2N.
    • Solve the system for c=1,f=0, a=0,b=20. alpha=beta=1 (cf. the problem above) i.e. compute the discrete solution, plot the graph. Compare with the shooting method i.e. the previous problem.
  • Solve the same problem -u''+c*u=f but with changed the right boundary conditions u(a)=alpha,u'(b)=beta - this condition for simplicity approximate by a backward difference (of order one unfortunately). Try different approach, namely add a ghostpoint x(N+1)=b+h and approximate the right bnd with the central diffrence i.e. uh(N+1)-uh(N-1)/(2h)=beta - compare the order of convergence... (naturally we have to assume that the differential equation is satisfied at x=b i.e. we have an extra equation: (-uh(N-1)+2uh(N)-uh(N+1)/(h*h)=f(b)).
  
  lab12.m - two examples of application of shooting method for solving BVPs - linear example u''=u u(0)=u(T)=1 - fails for large Tl and a nonlinear one u''=(3/2)u^2 u(0)=4;u(1)=1 - two solutions; an application of FDM for the first linear problem u''=u u(0)=u(T)=1 works for any T (reasonable)
• Simple control problem for IVP or BVP
  • A very simple control problem for linear elliptic BVP (heated rod - ery close to the simplest version of the project )
    Find s0 : F(s)=\int_b^1|y(s)-H(t)|^2 dt is minimal for s0 where y solves -y''=1 y(0)=s0;y'(1)=0; H- given function on [b,1] (0 < b < 1), (simplified project i.e. 1d instead 2d) Solve the problem by FDM method on the equidistant mesh x(k)=k*h k=0,..,N with h=1/N.
    Use backward difference to approximate the right bnd condition. This gives us a 1st order scheme.
    In order to get a 2nd order scheme add a ghostpoint x_{n+1)=1+h; and approximate the bnd cnd by the central difference u(N+1)-u(N-1)/(2h)=0 adding also the FDM approximation of the Diff equation at x=1 i.e. (-u(n-1)+2u(n)-u(n+1))/(h*h)=1 assuming that -u''(1)=1. (We have an extra value at the ghost point so we need an extra equation).
  • Find s0 : f(s)=\int_c^b|y(t;s)-h(t)|^2 dt is minimal for s0 where y solves y''=F(t,y,y'); y(0)=0;y'(0)=s; h- given function on [c,b] (a < c < b), F given e.g. F(x)=-sin(x) penduulum equation; One can replace IVP with BVP e.g. demanding that y solves y''=F(t,y,y'); y(0)=0;y'(b)=s; etc
  
  lab13.m - a solution of the control problem find s minimizing F(s)=\int_0.5^1 |u-H|^2 dx with u solving -u''=f u(0)=s u'(1)=0 - we solve the problem with f=-2 and H=Const=20; try different f and H. Requires fdmbvpsolmix.m script, see below.
  fdmbvpsolmix.m -FDM solver of the problem -u''=f u(a)=al u'(b)=be
• LAB: Section 6 in the lecture notes and a parabolic control problem
  • Section 6 in lecture notes: Charakterystyka pracy transformatora.(in Polish) we have ODE: u''+(w0^2)*u+b*u^3=(w/N)E cos(wt)
    N - no of coils of transformer, AC of frequency w/(2*\pi) and E- voltage, ; b,w0- parameters of a tranformer
    Plot the solution for E=165, w=120\pi, N=600, w0^2=83,b=0.14 on K=100 points, then check if solver is OK decreasing the tolerances (multiply them by 1e-4), if it is Ok then check the graph plotting for K=300,600. How to pick the right value of K? (1000,3000,5000?) And how to plot reasonable part of the grap e.g. 20,10,5% of the graph?
  • Parabolic PDE in 1D u_t=u_{xx} (0,T] x (0,L) - space 1D -rgt bnd : homogeneous Neumann u'(t,L)=0; left bnd cond : Dirichlet - u(t,0)=s(t), u(0,x)=u0(x) given. Solve by FDM using lsode() after space FDM dicretization or solve ODE by using implicit Euler scheme.
    Semidiscrete ODE system u(t,x(k))'- (u(t,x(k-1))+2(u(t,x(k))- (u(t,x(k+1)))/(h^2)=f(t,x(k)) k=1,..,N-1 u(t,x(0))=s(t); [u(t,x(N))-u(t,x(N-1))]/h=0; u(0,x(k))=u0(x(k)) k=1,..,N-1; x(k)=k*h h=L/N;
  • Simple control problem for time dependent PDEs (1D in space) Find s0(t) such that f(s)=\int_0^T \int_c^b|y(t,x;s)-H(t,x)|^2 dx is minimal for s0 where y solves y_t-y_{xx}=0 y(0,x)=y0(x) given, y_x(t,1)=0 and y(t,0)=s0(t) (control), H- given function on [c,b] (a < c < b), One can replace Dirichlet bnd by Robin/Neumann one for s0. - in practice we can look for s which is polynomial (constant, linear, quadratic, cubic etc) or a linear spline defined by values at discrete times at which we compute the values of u using lsode() - but those nodes of this linear spline should be independent of h (naturally we may assume that s a continuous function approximated by our scheme and look for all s(n*tau) - tau discrete time step) Take T=2, a=-1, b=1, c=0, H=20 (constant wit hrespect to all variables). Try first constant s0 (a constant polynomial) then s0 linear (the simplest is s0(t)=a+ ((b-a)/T)t -a is starting temeprature linearly changed to b (for time T). In general, for a polynomial s0 present the polynomial s0(t) by its values at equidistant points on [0,T] - this values may represent temperature values at those time points.
  • The same control problem but with Robin bnd i.e. -u'(t,a)+ u(a,b)=s(t) - u' can be approximated by the forward difference

Octave scripts, m-files,simple C files with e.g. examples of an use of BLAS/LAPACK

Projects

1. Optimal control (a stationary version in 1D discretized by FDM is easy and recommended by all of You) There is one project in two versions, i.e. optimal control of a heated room, described in details (there 2 versions: stationary and unstationary with many possible space discretizations method - I would describe 2 in details during lectures). I have other projects - ask for details).
2. (medium difficult) ODE model of golf putting - in the paper a model of golf putting on the green is presented. The goal of project using pctave ODE solver (e.g. lsode() or oder45()) solve the model for examples of green surface presented in the paper - repeat the computations - compare with the tables from the paper.
   In the optimal control project I simplified necessary tests (it is enough to test the linear case - cf. describtion for details)

   ---

   
   Return to my homepage
   
   

   ---

   Last update : May 28, 2019
   
   

   ---