Root Finding Bracketing Methods 1. Bisection Method - Isolate all x (f(x)=dasdasd) - Plot f(x).

Find an interval [a,b] where 1 root exists. Interval such that f is monotonic in [a,b] and f(a)*f(b)<0 - Find midpoint of [a,b] and label it as c. Evaluate f(a), f(b) and f(c) - If f(a)*f(c)<0, then set new bounds to [a,c]. If f(b)*f(c)<0, then set new bounds to [b,c] -Iterate depending on tolerance. Let say tolerance = 0.01. The smaller the tolerance, the more accurate. TOLERANCE > range of bounds. Usually 0.0001. -a,b,c can be a root. Usually c. Instead of having a tolerance, you can also set the maximum # of iterations. format long %bisection f = inline('2^x-x^2-2'); a = input('a: '); b = input('b: '); tol = 0.0000001; while tol <= abs(f(a)-f(b)) c =0.5*(a+b) if f(c)*f(a) < 0 b=c; end if f(c)*f(b) < 0 a=c; end end disp(c); 2. Regula-Falsi Method (False Position) - Isolate all x (f(x)=dasdasd) - Plot f(x). Find an interval [a,b] where 1 root exists. Interval such that f is monotonic in [a,b] and f(a)*f(b)<0 - Solve c using c = b- f(b)(b-a)/(f(b)-f(a))=b-f(b)/slope Derivation from similar triangles: f(b)/f(b)-f(a)=(b-c)/(b-a) -- If f(a)*f(c)<0, then set new bounds to [a,c]. If f(b)*f(c)<0, then set new bounds to [b,c] -Iterate depending on tolerance. Let say tolerance = 0.01. The smaller the tolerance, the more accurate. TOLERANCE > range of bounds. Usually 0.0001. -a,b,c can be a root. Usually c. -faster than Bisection format long %bisection f = inline('2^x-x^2-2'); a = input('a: '); b = input('b: '); tol = 1000; while tol >= 0 slope=(f(b)-f(a))/(b-a); c = b - f(b)/slope if f(c)*f(a) < 0 b=c; end if f(c)*f(b) < 0

a=c; end tol = tol-1; end disp(c); format long %bisection f = inline('2^x-x^2-2'); a = input('a: '); b = input('b: '); tol = 0.00001; d=[abs(a-b)]; while tol <= min(d) slope=(f(b)-f(a))/(b-a); c = b - f(b)/slope; d=[abs(a-b), abs(b-c), abs(a-c)]; if f(c)*f(a) < 0 b=c; end if f(c)*f(b) < 0 a=c; end end disp(c); Problems in Bracketing -Knowing where to start -Slow converges -Problem for small non-monotonic -Not good for eqn's for more than one root. Root Finding Open Methods 1. Fixed-Point Iteration -Isolate x on the left side. Declare whole right side as g(x). (form x = g(x)) -Greatly dependent on how x is isolated. (e.g. look for undefined cases!!!) -Make an initial guess, then iterate (g(g(g(g...(x))) until converges -Converges if |g'(xo)|=1 format long %fixedpointtrululu %'2^x-x^2-2' f = inline('log2(x^2+2)'); x=input('x: '); tol = 0.00000001 while tol <= abs(f(x)-x); x=f(x) end disp(x); 2. Newton-Raphson Method -Isolate all x -Plot f(x). Make an initial guess close to the root -x(i+1)-f(x(i))/f'(x(i)) -Tolerance>delta iterations Problem: Might diverge base on initial guess, not readily differentiable, sometimes guess is on relative extrema (diff = 0)

format long %newtonrhapson %'2^x-x^2-2' f = inline('2^x-x^2-2'); df = inline('2^x*log(2)-2*x'); x=input('x: '); tol = 0.00000001; y=0; while tol <= abs(x-y) y=x x=x-(f(x)/df(x)) end disp(x); Problem 1: Multiple Roots clear,clc p = 'log(x^2+1)-exp(0.4*x)*cos(pi*x)'; f = inline(p); tol = 0.000001; a = input('a:'); b = 0; c = (a+b)/2; while tol <= abs(f(a)) c =0.5*(a+b); if f(c)*f(a) < 0 b=c; end if f(c)*f(b) < 0 a=c; end end disp(c); Problem 2: Crossed Ladders clear,clc h=input('x: '); k=input('y: '); c=input('c: '); f = inline('1/sqrt(h^2-x^2)+1/sqrt(k^2-x^2)-1/c'); df = inline('x*(1/(h^2-x^2)^1.5+1/(k^2-x^2)^1.5)'); x=input('guess: '); y=0; tol = 0.0000000001; while tol <= abs(y-x) y=x; x=x-f(c,h,k,x)/df(h,k,x); end disp(x); Problem 3: Mechanics of Deformable Bodies clear,clc L = 600; k = 2.5/(120*50000*30000*L); p = '-5*x^4+6*L^2*x^2-L^4'; q = 'k*(-x^5+2*L^2*x^3-L^4*x)'; f = inline(p); g = inline(q); tol = 0.000001; a = input('a:');

b = input('b:'); c = (a+b)/2; while tol <= abs(a-b) c =0.5*(a+b); if f(L,a) == 0 || f(L,b) == 0 || f(L,c) == 0 break; end if f(L,c)*f(L,a) < 0 b=c; end if f(L,c)*f(L,b) < 0 a=c; end end %disp(['Point of max deflection: ', num2str(c)]); %disp(['Maximum deflection: ', num2str(g(L,k,c))]); fprintf('Point of max deflection: %.3f\n', c); fprintf('Maximum deflection: %.3f\n', g(L,k,c)); Numerical Integration a. Trapezoidal Rule b. Simpson's Rule -Recall the principle of Riemann sum :) -Some functions are difficult to integrate Trapezoidal Rule - Recall area of trapezoid: I = (b-a)(f(a)+f(b))/2 - Approximation of trapezoids, to improve the accuracy, we increase the number of trapezoids/segments: h= (b-a)/n = step size, n is now number of trapezoids I is now = h(f(ix)+f(ix+1))/2 AREA = summation of areas format long f=inline('abs(x)'); x=-4; y=4; l=abs(x-y); n=100; integral=0; h=l/n; for i = 0:n-1 integral=h/2*(f(x+i*h)+f(x+(i+1)*h))+integral; disp(integral); end Simpson's Rules - Based on Newton-Cotes quadratures Simpson's 1/3 Rule - Parabolic approximation: Area of an Approximate Parabola I = (b-a)(f(xo)+4(x1)+f(x2))/6 - Uses three points, thus it uses two portions. Number of steps must be EVEN. h= (b-a)/n I = b-a (f(xo)+4sum(i=1,3,5 to n-1 of f(x1)) + 2sum(i=2,4,6 to n-2 of f(x2))/6 format long f=inline('1-exp(-2*x)'); x=0;

y=4; l=abs(x-y); n=10; integral=0; h=l/n; for i = 0:2:n-2 integral=h/3*(f(x+(i)*h)+4*f(x+(i+1)*h)+f(x+(i+2)*h))+integ ral; disp(integral); end format long f=inline('1-exp(-2*x)'); x=0; y=4; l=abs(x-y); n=10; integral=0; b=linspace(x,y,n+1); h=l/n; for i = 1:2:n-1 integral=h/3*(f(b(i))+4*f(b(i+1))+1*f(b(i+2)))+integral; disp(integral); end

end a=input('x: '); format long f=inline(['exp(-t)*t^(',num2str(a),'-1)']); x=0; y=100; l=abs(x-y); n=3300; integral=0; b=linspace(x,y,n+1); h=l/n; for i = 1:3:n-2 integral=3*h/8*(f(b(i))+3*f(b(i+1))+3*f(b(i+2))+f(b(i+3)))+i ntegral; end MACHINE EXERCISES Problem 15: Root Finding and Integration clear,clc format long df = inline('1/sqrt(2*pi)*exp(-t^2/2)'); N=input('N: ');%area g=0.5;%guess tol = 0.00000001; y=0;%previousvalue %integrationconditions x=0; n=100; while tol <= abs(g-y) y=g; l=abs(x-g); h=l/n; integral=0; for i = 0:n-1 integral=h/2*(df(x+i*h)+df(x+(i+1)*h))+integral; end g=g-((integral-N)/df(g)); end fprintf('%.4f\n', g); Sample Machine Problem 3 (You can use matrices in integration. ITS FASTER.) clear,clc format Nvalue=2:20 for k=1:length(Nvalue) N=Nvalue(k); f=inline(['t.^',num2str(N),]); x=0; y=1; l=abs(x-y); n=5000000; h=l/(n-1); fmatrix=f(x:h:y); intmatrix=zeros(1,n);

Simpson's 3/8 Rule Area of an Approximate Parabola: I = (b-a)/8(fo+3f1+3f(2)+1)/8 - n must be divisible by THREE format long f=inline('1-exp(-2*x)'); x=0; y=4; l=abs(x-y); n=33; integral=0; h=l/n; for i = 0:3:n-3 disp([x+(i)*h,x+(i+1)*h,x+(i+2)*h,x+(i+3)*h]) integral=3*h/8*(f(x+(i)*h)+3*f(x+(i+1)*h)+3*f(x+(i+2)*h)+f( x+(i+3)*h))+integral; disp(integral); end format long f=inline('1-exp(-2*x)'); x=0; y=4; l=abs(x-y); n=333; integral=0; b=linspace(x,y,n+1); h=l/n; for i = 1:3:n-2 integral=3*h/8*(f(b(i))+3*f(b(i+1))+3*f(b(i+2))+f(b(i+3)))+i ntegral; disp(integral);

integral=0; for i = 0:n-2 intmatrix(i+1)=h/2*(fmatrix(i+1)+fmatrix(i+2)); end intmatrix=1./sqrt(fliplr(cumsum(fliplr(intmatrix)))); for j = 0:n-3 integral=h/2*(intmatrix(j+1)+intmatrix(j+2))+integral; end integral=integral/sqrt(pi); test=sqrt(N+1)*gamma((N+2)/(N+1))/gamma((N+3)/(2*N +2)); disp(N); fprintf('%.2f\n',integral); disp(test); end Numerical Differentiation 1. Forward Difference 2. Backward Difference 3. Central Difference - Using small intervals between values of x, the derivative can be approximated using limits: Given the interval [a,b], and all values xn is in the interval: Forward Difference: f'(xi) = (y(i+1)-yi)/(x(i+1)-xi) Applicable at [a,b) Backward Difference: f'(xi) = (y(i)-y(i-1))/(x(i)-x(i-1)) Applicable at (a,b] Central Difference: f'(xi)= (y(i+1)-y(i-1))/(x(i+1)-x(i-1)) Applicable at (a,b) MACHINE EXERCISE

if i > 1 trapezoid(i-1)=0.4*(deriv(i)+deriv(i-1))*0.5; end end t=sum(trapezoid)/(V*60/1000) + (N-2)*3; disp(t); Various Problems: FLAMES clear,clc %flames flames = 'FLAMES'; A=input('name 1: ', 's'); B=input('name 2: ', 's'); c=0; d=intersect(A,B); A=[A,B]; if isempty(d) disp('NO'); else for i=1:length(d) c=length(find(A == d(i)))+c; end counter=mod(c,6); if counter == 0; counter=6; end disp(flames(counter)); end Frankline Temperature x1 = input('x1: '); y1 = input('y1: '); x2 = input('x2: '); y2 = input('y1: '); T = input('T: '); A = (y1-y2)/(x1-x2); B = y1 - x1*A; frankline = A*T+B; fprintf('%.4f\n', frankline); Residue Function - Long Version

Problem 17: MRT N = input('N: '); V = input('V: '); x = 0:0.4:10; y = [0.200,0.049,0.015,0.052,0.275... 0.670,1.241,1.982,2.886,3.938... 5.122,6.420,7.810,9.275,10.794... 12.349,13.921,15.494,17.052,18.581... 20.069,21.504,22.878,24.181,25.408,26.553]; deriv = zeros(1,length(x)); trapezoid = zeros(1,length(x)-1); for i = 1:length(x) if i == 1 deriv(i) = (y(i+1)-y(i))/(x(i+1)-x(i)); elseif i == length(x) deriv(i) = (y(i)-y(i-1))/(x(i)-x(i-1)); else deriv(i) = (y(i+1)-y(i-1))/(x(i+1)-x(i-1)); end deriv(i)=sqrt(1+deriv(i)^2); clear, clc; p = input('P: '); q = input('Q: '); r = input('R: '); if p == q && p == r A = 0; B = 0; C = 1; elseif p == q || q == r || p == r s = unique([p,q,r]); B = 1/(s(1)-s(2)); C = 1/(s(1)-s(2))^2; A = -C; else A = 1/(p-q)/(p-r); B = 1/(q-r)/(q-p); C = 1/(r-p)/(r-q); end fprintf('%.4f ', [A, B, C]);

Curve-Fitting Techniques Linear Regression - Sum of squares of residuals (SSr) must be LEAST: - R^2 value must be closest to ONE: SSr = Summation(Y Yfit)^2, where Yfit is the Y value in the linear fit. SSt = Summation(Y-Ymean)^2, where Ymean is the mean of actual Y values R^2 = 1 (SSr/SSt), if SSr is least, then R^2 value is greatest, as SSr approaches zero, R^2 approaches 1. Linear models for non-linear equations: 1. Exponential: lny= lnb + mx 2. Power: logy = logb + mlogx 3. Saturation-Growth Rate: 1/y = b/a(1/x) + 1/a MATLAB has a built-in function that fits a least-squares nth-order polynomial to data p = polyfit(x,y,n) p = vector containing coefficients x and y = vector containing data n = desired polynomial order You can use the p vector to evaluate y-fit values: y = polyval(p,x) y = yfit, p = vector of coefficients, x = x-data MATLAB follows a specific format for a polynomial: [p1, p2] = p1x + p2 [p1, 0, p2] = p1x^2 + p2 [p1,p2,p3,0,p4,2,p5] = p1x^6 + p2x^5 + p3x^4 + p4^2 + 2x + p5 General Linear form: y = a0f0(x) + a1f1(x) + + amfm(x) where a = coefficients, f(x) = sub-functions Example: y=1+x a0 = 1, f0 = 1, a1 = 1, f1 = x y = 3 + 2x + x^3 a0 = 3, f0 = 1, a1 = 2, f1 = x, a2 = 0, f2 = x^2, a3 = 1, f3 = x^3 Procedure: 1. Enter the data to be fit (y and x) as column vectors 2. Select a model you wish to fit your data. 3. Create your z matrix [f0(x0). f1(x0),;f0(xn), f1(xn)fn(xn)]; 4. Solve the following: A = Z*Z; B = Z*y;

5. Solve for matrix containing coefficients: X = A\B; (inverse matrix), X is column vector Example: Find a quadratic equation: if the x-values are 0-5, then: f0 = 1, f1 = x, f2 = x^2 z = [1, 0, 0^2;1, 1, 1^2; 1, 2, 2^2;1, 3, 3^2; 1, 4, 4^2; 1, 5, 5^2;] z = [1, 0, 0; 1,1,1; 1,2,4; 1,3,9; 1,4,16; 1,5,25;] Given an equation y = mx+b, the equations needed to solve the coefficients are: m = N(sum(xi*yi))-sum(xi)*sum(yi)/(N(sum(xi^2))(sum(xi))^2)) b = (sum(yi)-m(sum(xi)))/N N = number of data points xi and yi = x-data and y-data given in ith point Problem 24: Michaelis-Menten Kinetics x=1./(0.5:0.5:5)'; y=1./[4.023,4.121,4.152,4.158,4.163,4.179,4.182,4.184, 4.188,4.191]'; %f1(x)=1 %f2(x)=x Z=[ones(length(x),1), x]; A=Z'*Z; B=Z'*y; solution=A\B; Vmax=1/solution(1); Km=solution(2)*Vmax; fprintf('%5d \n', Vmax); fprintf('%5d \n', Km); disp(Km); Problem 23: Gompertz Tumor Growth Model This Problem uses lsqcurvefit function. function F = tumor(x, xdata) F = x(1)*exp(log(580.23/x(1))*exp(-x(2)*xdata)); end format shortG; clear,clc xdata=0:2:30; ydata=[580.23, 865.34, 1170.25... 1529.30, 1899.68, 2273.95... 2623.62, 2938.75, 3246.85... 3510.31, 3740.78, 3947.47... 4118.59, 4252.85, 4382.22... 4489.60]; x0=[100;10];%startingguess [x, resnorm]=lsqcurvefit(@tumor, x0, xdata, ydata); i1 = input('x1: '); %i2 = input('x2: '); xdata = i1; ydata = tumor(x, xdata); disp([xdata, round(ydata)]);