Engineering Computation Workshop 9 === ###### tags: `comp20005` `workshop` `c` starting at 6:20 --- ### Assignment 1 Feedback (General) - Function > 60 lines --- - Function takes more than 3 arguments --- - No Comments --- - Commenting on (if (x > 0)) // check if x is greater than 0 --- - Writing the whole thing in main --- - Writing anything in main (except func calls) --- - Super long variable and function names --- - 80 char line limit --- - duplicate code (printing or across functions) --- ### Root Finding --- ### Binary Search How do we find an element in a sorted sequence in $O(logn)$ as opposed to linear scan $O(n)$? --- ### Binary Search By throwing away half of the options every check! --- ### Binary Search  --- ### Binary Search General algo ``` A = sorted_array, x = element_to_search l = left_index, r = right_index function bsearch(A, x, l, r): while (l <= r): m = middle(l,r) if (A[m] == x): return m if (A[m] < x): l = m + 1 else: r = m - 1 return -1 ``` --- ### Binary Search ```cpp int bsearch(int* A, int x, int l, int r) { if (r < l) return -1; int m = l + (r - l) / 2; // overflow safe mid if (A[m] == x) // found it return m; if (A[m] > x) // its on the left side return bsearch(A, x, l, m-1); return bsearch(A, x, m+1, r); // right side } ``` --- ### Bisection Root Finding Can we apply the same logic to finding the result of a function? --- ### Bisection Root Finding Yes!  --- ### Bisection Root Finding General Function Binary Search on monotonic integer functions ```cpp #define NONE -1 #define INC 1 int f(int x); int bsearch(int x, int l, int r){ if (r > l) return NONE int m = l + (r - l) / 2; if (m == x) return m; if (f(m) > x) return(x, l, m-INC) return(x,m+INC,r) } ``` --- ### Bisection Root findin When we want to find roots? ```cpp #define EPS 10e-4 #define LIMIT 10e5 int d(int x); int bisect(double x, double l, double r, int c){ double m = l + (r - l) / 2; if (c == LIMIT || l > r) return m; if (d(m) * d(l) < 0) //root to left return(x, l, m); return(x,m,r); // root to right } bisect(0, 10, -10, 0) ``` --- ### Newton Raphson Using the slope of the function, we can compute an appox. root using the tangent of the curve $$x_{i+1} = NR_{x_i} = x_i - \frac{f(x_i)}{f'(x_i)}$$ --- ### Newton Rapson ```cpp #define EPS 10e-4 double f(double x); // func double df(double x); // deriv func double NR(){ int i, max_i; long double h, x0, x1, allerr; for (i = 1; i < max_i; ++i){ h=f(x0)/df(x0); x1=x0-h; if (fabs(h) < EPS) return x1; x0 = x1; } return -1; } ``` --- ### Integration --- ### Trapezoid  --- ### Trapezoid Estimate by breaking the interval into n stpes, of size $h=\frac{(x_2 - x-1)}{n}$. Then we can compute the area as: --- $\sum_{i=0}^{n-1}{h\frac{f(x_1 + h\cdot i) + f(x_1 + h \cdot (i+1))}{n}}$ $=\frac{h}{2}(f(x_1) + 2 \cdot \sum_{i=1}^{n-1}{f(x_1 + h \cdot i) + f(x_1 + h \cdot n)})$ --- ### Trapezoid ```cpp double f(double x) double trp(double x1, double x2, double n) { double h = (x2-x1)/n; double s = f(x1) + f(x2); for (int i = 1; i < n; ++i) s += 2 * f(x1+i*h); return (h/2)*s; } ``` --- ### Simpsons Fitting quadratic to triples of points. Given $x_1,x_2,x_3$ are such that $x_2-x_3=x_3-x_2$ and $P(x) = p\cdot x^2 + q\cdot x + r$ is fitted through $(x_1,f(x_1)), (x_2,f(x_2), (x_3,f(x_3))$ then we get an easily solvable method. --- $\frac{h}{3} (f(x_0) + f(x_0+h\cdot n)$ $+4\cdot \sum_{i=1,3,..}^{n-1}{f(x_0 + h \cdot i)}$ $+2\cdot\sum_{i=2,4,..}^{n-2}{f(x_0+h\cdot i)})$ --- ### Simpsons ```cpp double simp(double a, double b, int n) { double h = (b-a)/n; double res = 0; for (int i=0;i<=n;++i){ if (i == 0 || i == n) res += f(a + i * h); else if (i % 2 == 0) res += 4 * f(a + i * h); else res += 2 * f(a + i * h); } return res * (h/3); } ``` --- ### Interpolation No function, given points and have to estimate roots / integral - Linear - Newton --- ### Linear  --- ### Linear Piecewise linear interpolation $\hat{f_1}(x) = f(x_i) + (x-x_i)\frac{f(x_i+1)-f(x_i)}{x_i+1 - x}$ --- ### Newton Fitting a lesser degree (by 1) polynomial through n points --- ### Newton $\hat{f_2}(x) = b_1 + b_2(x_i) + b_3(x-x_1)(x-x_2)$ --- ### Newton Fitting to $px^2 + qx + r$ we have $p = b_3$ $q = b_2-b_3x_1 - b_3x_2$ $r=b_1-b_2x_1+b_3x_1x_2$ --- ### Systems of Linear equations - Guassian - LU - LUP --- ### Guassian ``` Input: For N unknowns, input is an augmented matrix of size N x (N+1). One extra column is for Right Hand Side (RHS) mat[N][N+1] = {{3.0, 2.0,-4.0, 3.0}, {2.0, 3.0, 3.0, 15.0}, {5.0, -3, 1.0, 14.0} }; Output: Solution to equations is: 3.000000 1.000000 2.000000 Explanation: Given matrix represents following equations 3.0X1 + 2.0X2 - 4.0X3 = 3.0 2.0X1 + 3.0X2 + 3.0X3 = 15.0 5.0X1 - 3.0X2 + X3 = 14.0 There is a unique solution for given equations, solutions is, X1 = 3.0, X2 = 1.0, X3 = 2.0, ``` --- ### Guassian 1. Find the the kth pivot by swaping rows, to move the entry with the largest abs value to thee pivot position 2. For each row below the pivot, calc the factor f which makes the kth entry zero, for every element in the row subtract the fth multiple of the corresponding element in the kth row 3. Repeat for each unknown --- ### Guassian Elimination $3a + 2b + 4c = 3$ $2a + 3b + 3c = 15$ $5a - 3b + 1c = 14$ Matrix form: $\begin{matrix}[ 3 & 2 & -4 & 3 ]\\ [2 & 3 & 3 & 15 ]\\ [5 & -3 & 1 & 14 ] \end{matrix}$ --- ### Guassian Elimination Swap R1,R3 $\begin{matrix}[5 & -3 & 1 & 14 ]\\ [2 & 3 & 3 & 15 ]\\ [ 3 & 2 & -4 & 3 ] \end{matrix}$ --- ### Guassian Elimination R2 = R2 - R1/2.5 $\begin{matrix}[5 & -3 & 1 & 14 ]\\ [0 & 4.2 & 2.6 & 9.4 ]\\ [0 & 3.8 & -4.6 & -5.4 ] \end{matrix}$ --- ### Gaussian Elimination Second Pivot $\begin{matrix}[5 & -3 & 1 & 14 ]\\ [0 & 4.2 & 2.6 & 9.4 ]\\ [0 & 0 & -6.9 & -13.8 ] \end{matrix}$ --- ### Gaussian Elimination Result $\begin{matrix}[5 & -3 & 1 & 14 & =3]\\ [0 & 4.2 & 2.6 & 9.4 & =1]\\ [0 & 0 & -6.9 & -13.8 & =2] \end{matrix}$ ---