Scientific Computing 2015/16

Konsultacje:

References

1. P. Krzyzanowski, Obliczenia inzynierskie i naukowe. Skuteczne, szybkie, efektowne, Wydawnictwo Naukowe PWN, 2011
2. P. Krzyzanowski, Obliczenia naukowe. Html lecture notes (in Polish). Pdf version is alos available there. University of Warsaw. 2010.

1. Lab 1 - 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 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)
2. Lab 2- scripts: example matbasic.m write on command line of octave'a: matbasic.
   1: using vander() and polyval() - 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().
   2: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 i 30 nodes.
   3- do same as in prob 2 but for Chebyshev points (linearly scaled roots of T_n(x)=cos (n*arc cos(x))).
   lab2.m - a few solutions
3. Lab 3 and 4 Solve nonlinear equations (fsolve and fzeros) for simple equations x^2-2=0 etc exp(x)+x=0 . Global variables can be useful.... Prob 1: Solve x*(1+0.5sin(x))=0; (x-1)^2=0; Prob 2: 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. Prob 3: plot the graph of the inverse of sin(x) on [0,pi] or x^ on [0,2] (compare with arcsin(x) or \sqrt(x)) etc or whe ndo not know the direct formula for the inverse f.: (x+1)*exp(x^2) on [1,3] or x+cos(x). Prob 4: 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]
   Serious project we want to plot the inverse of a given function - see formulas in Chapter 3 in lect. notes.
   nonlin16.m - some simple examples of solving nonlinear equations (fsolve and fzero)
4. Lab 5 and 6 : Numercal classics, sparse matrices, loops, integration etc : LU QR decompositions, LLSP, how to compute the inverse of a matrix, condition number, norms, SVD, null space, R(A) - how to get orhonormal basis of a system of vectors etc Sparse matrices global and static (persistent a=0) variables, .
   • Prob : compute a 3-diagonal matrix without forming a full format matrix and a sparse format matrix (parameters: a -diagonal,b-subdiagonal,c-superdiagonal i.e.vectors) m-files- how to define functions, how to write help to our functions, Compare time of creating full/sparse formats of 3diag matrix wit h2 on main diagonal and -1 on sub- superdiagonals, then compare times of solving linear systems with this matrix with both formats.
     on16sparse.m -some solutions for sparse 3diag matrices
   • nargin, nargout, varargin etc Subfunctions, functions within functions etc
   • Conditionals, loops etc if else, switch, pętle while, do until, logicals (&& and, || - or, !-not, != - not equal == equality), in particular x||y for x,y, vectors/matrices
   • Prob : implement without loops a hut 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(), ), całki z osobliwościami, <
   • 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 a 2d volume of an elipsis 2x^2+3y^2-1=0 (y(x) is given implicitely... )
5. Because of a computer science high school competition - last lab was cancelled - we can do 1-2 problems from the previous lab. 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 uop equation: y'=y*y; y(0)=1 on [0,1/2] or [0,3]. Linear penduulum equation: y''=-y z y(0)=1;y'(0)=0 non[0,2*pi](solution: cos(x)). 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 2 NPlot 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 3 Find s0 such that: y'=y*(1-y)*sin(y*y+x); y(0)=s0 satisfies y(1,s0)=0.5. Plot the graph of F(s)=y(1,s) and y(t,s0) for t in [0,1].
   • Prob 4 Solve boundary value problem: : u''+u*u=f(t) u(0)=u(1)=0 for different f e.g. f=0; 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 i
6. Earthquake model, ODEs cont. Solving ODEs cont.
   • Earthquake model. of a building Stiff floors - with elastic connections between them and the zero floor wit hthe ground, assume that the i+1 th floor of the mass m_{i+1} which moves relativly to ith floor: x_{i+1}-x_{i} then a restoration force k_i(x_{i+1}-x_i) (proportional to the relative shift) acts on the both i.e. on the i-th acts 2 forces: k_i(x_{i+1}-x_i) and -k_{i-1}(x_{i}-x_{i-1}), on the zero floor (ground floor) acts analogous restorative force proportional to the shift relative to the ground: i.e. -k_0x_1 and the restorative force between zero floor and 1st floor During than earthquake additionaly the force f(t)=Gsin(gamma t) acts on the ground floor (level 0) e.g .gamma=3 we have ODE z MX'=SX+F for M diagonal matrix n+1 x n+1 (masses of floors on the diagonal), S 3diagonal symmetric matrix such that on the main diagonal we have -k_i-k_{i-1} (k_{n+1}=0=0) sub- super digonal = (k_1,..,k_n) F=(f(t),0...,0)' z f_1(t)=Gsin(gamma*t) - G<>0 during earthquake (in practise G should of the form G(t)=G exp(-a|t-t_0|) a - positive large
     • Write ODE: NMX''=SX+ CX' + F(t). (F =(f(t),0..,0)^T, C - diagonal matrix - air resistance. Write functions computing matrices: A=M^{-1}S, M^{-1}C and 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 witout 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}....
7. Simple control problem for IVP or BVP
   • 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
   • Another very simple control problem for linear elliptic BVP (heated rod)
     Find s0 : f(s)=\int_c^b|y(t;s)-h(t)|^2 dt is minimal for s0 where y solves -y''=f y(0)=s0;y(1)=0; h- given function on [c,b] (a < c < b), (simplified i.e. 1d instead 2d - project) Sove the problem by FDM method.
   • 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?
   ode16-1.m - some solutions for odes
   earthq16.m - some solutions earthquake model
   cp1d.m - a script with the solution of 1d control problem -u''=f u(0)=s u'(2)=0 we want to find s0 such that u(x;s0) minimizes : \int_1^2 |h(x)-u(x;s)|^2 dx 1/ test for f=0 (then u(x;s)=s) an h=0.5 the naturally s0=0.5
   How to organize scripts i.e. path,addpath etc.
8. Control problem - 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) - control in pracitise let s be a lineatr spline (lpiecewise linear between times; or a polynomial of set degree e.quadratic. u(0,x)=u0(x) given.
   • Finished previous problems (again it would take 2 labs so we can only try if time allows) Solving PDEs cont. Control problem - parabolic PDEs
   • 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.
     cpp1d.m - a script with the solution of 1d control problem (unfinished) u_t-u_{xx}=0 u(t,0)=s(t) u'(t,2)=0 we want to find s0(t) such that u(x;s0) minimizes : F(s)=\int_0^T \int_1^2 |h(t,x)-u(t,x;s)|^2 dx dt - h is a given function e.g. h constant - - in practise 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()
     in the script we have a function solving parabolic problem for any s(), it remains to write the objective function i.e. approximation of F(s) and to minimize it finding s0
   • The same control problem but with Robin bnd i.e. -u'+ u=s(t) approximated by FDM (u_0(t)-u_1(t))/h+u_0(t)=s(t) --> u_0=(s+u_1)/(1+h) for all t
9. C language and using fortran libs i.e. e.g. BLAS in C/C++
   1. Prob 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
   2. Prob 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].
      vecmatsave.c - a simplec C code writing vectors/matrices to text files
      Homework: 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.
   3. Prob Alokacja pamięci. 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 : apointer 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
   4. Prob 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).
   5. 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, adds two matrices (in that format), multiplies two matrices, scales a matrix, creates matrix of given dimension (allocates memory), destroys matrix (frees memory) etc
   6. 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.
   7. 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.g 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 lmath lib)
      tblas.c - a simple C code using BLAS dnrm2() function - and dynamically allocating memory (using malloc(); free() functions)
10. Using fortran libs cont.: matrix operations in BLAS, solving linear systems and overdetermined linear systems (LLSP) in LAPACK in C/C++
    1. Prob 1 Write 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. Using BLAS DGEMV(): dgemv.f compute Ax for A=[1 2;-1 1] and x= [0;1], [1;0], [1;1].
    2. Prob 2: Solve system Ax=f for A=[ 2 -2; 1 1] and f=[0; 2] (rozw [ 1;1]) using LAPACK function DGESV(): dgesv.f
    3. Prob 3: 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.
    4. Prob 4: Find LAPACk function solving the tridiagonal system and use for the 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=100.
    5. Prob 5: 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]).

---

A few octave scripts, m-files, source C files, LAPACK functions

Link to Octave
octave-forge - octave extensions


Octave scripts, m-files

matbasic.m - a script with basic octave commands we can edit it : gedit matbasic.m
lab2.m - a few solutions of lab2
nonlin16.m - some simple examples of solving nonlinear equations (fsolve and fzero)
on16sparse.m -some solutions for sparse 3diag matrices
sesja4.m - a few simple examples on solving nonlinear equations and integration in octave
ode16-1.m -some solutions for odes (simple odes in octave 1d or 2d, blow-up real and false, stiff solved by explicit or implict methods in lsode)
sesja5.m - solving ODEs - 2 simple examples


---

Wyklad - Środowko programisty gdzie można znależć materiały o języku C jak i Makefile'ach (in Polish)

---

Project (in Polish). There is one project, 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).
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 : August 30, 2016



---