Here’s a brief MATLAB class containing methods for performing numerical analysis. Also available on GitHub.
General Usage
From a file in the same directory as NumUtils.m call:
NumUtils.MethodName(args)
Here is NumUtils.m:
classdef NumUtils % Defines numerical analysis utility functions. % % Methods: % AB2 ------------- Adams-Bashforth 2-step (AB2) LMM for solving IVP. % Bisection ------- Bisection Method. % CentralDiff ----- Central Difference for approximating f'(x0). % EstimateFPIter -- Estimates the number of iterations required for % convergence of Fixed Point Iteration algorithm. % EulersMethod ---- Euler's Method for solving an IVP. % FivePointMidpoint - Five-Point Midpoint Formula for approximating % f'(x0). % FixedPointIter -- Fixed Point Iteration. % ForwardDiff ----- Forward Difference for approximating f'(x0). % GaussSeidelMethod - Gauss-Seidel method for iteratively solving a % linear system of equations. % GramSchmidt ----- Construct orthogonal polynomials w/ Gram-Schmidt. % Jacobian -------- Symbolically calculates Jacobian for a system. % JacobisMethod --- Jacobi's method for iteratively solving a % linear system of equations % Lagrange -------- Generate Lagrange interpolating polynomial. % LLS ------------- Constructs linear least squares poly. coeffs. % LogB ------------ Calculates log(X) with base B. % NaturalCubicSpline - Calculates the natural cubic spline for f. % NewtonsMethod --- Newton's method for root finding problem % NewtonsMethodForSystems - Newton's Method for iteratively solving a % nonlinear system of equations F(x)=0. % QuasiNewton ---- Quasi-Newton Method for root finding problem % using global Bisection Method and local Newton's. % QuasiSecant ---- Quasi-Secant Method for root finding problem % using global Bisection Method and local Secant. % RK2 ------------- Runge-Kutta 2-step (RK2) for solving IVP. % SecantMethod ---- Secant method for root finding problem % TaylorPoly ------ Symbolically calculate the first N terms of a % Taylor Polynomial. % TaylorPolyNTerm - Symbolically calculate the Nth term of a Taylor % Polynomial. % ThreePointEndpoint - Three-Point Endpoint Formula for approximating % f'(x0). % ThreePointMidpoint - Three-Point Midpoint Formula for approximating % f'(x0). % TruncationError - Symbolically calculate the truncation error % associated with Taylor Polynomial approximation. % TruncationErrorLagrange - Symbolically calculate the truncation % error associated with Lagrange % Interpolating Polynomial approximation. % % Usage: NumUtils.FunctionName(args) methods(Static) %% AB2 function [t_vec,y_vec] = AB2(f,a,b,alpha,h,onestep) % Uses the Adams-Bashforth 2-step (AB2) linear multistep method % for solving an IVP. % % Input: f = (y'(t) = f(t,y)) % a = Start point of time interval. % b = End point of time interval. % alpha = y(a) % h = Timestep % onestep = One-step method to use for computing y1. % Select forward_euler, backward_euler, or rk2. % % Output: t_vec = Column vector of mesh points (each timestep) % y_vec = Approximated solution at each mesh point % Initialize output vectors t_vec=a:h:b; y_vec=zeros(length(t_vec),1); y_vec(1)=alpha; % Use onestep method to get y_1 switch onestep case 'forward_euler' [~,w]=NumUtils.EulersMethod(f,a,a+h,h,alpha,'forward'); % Store result of y_1 y_vec(2)=w(2); % Iteratively calculate solution at each mesh point for ii=3:length(t_vec) f_term=-f(t_vec(ii-2),y_vec(ii-2))+3*f(t_vec(ii-1),y_vec(ii-1)); y_vec(ii)=y_vec(ii-1)+(h./2).*f_term; end case 'backward_euler' [~,w]=NumUtils.EulersMethod(f,a,a+h,h,alpha,'backward'); % Store result of y_1 y_vec(2)=w(2); % Iteratively calculate solution at each mesh point for ii=3:length(t_vec) f_term=-f(t_vec(ii-2),y_vec(ii-2))+3*f(t_vec(ii-1),y_vec(ii-1)); y_vec(ii)=y_vec(ii-1)+(h./2).*f_term; end case 'rk2' [w]= NumUtils.RK2(f,[a;a+h],alpha); % Store result of y_1 y_vec(2)=w(2); % Iteratively calculate solution at each mesh point for ii=3:length(t_vec) f_term=-f([t_vec(ii-2),y_vec(ii-2)])+3*f([t_vec(ii-1),y_vec(ii-1)]); y_vec(ii)=y_vec(ii-1)+(h./2).*f_term; end otherwise disp('Error selecting onestep method. Please use ' ... + 'forward_euler, backward_euler, or rk2.') end end %% Bisection function [p,iter,relerr,p_all,iter_all,relerr_all] = Bisection(f,a,b,tol,maxiter,verbose) % Function implementing the bisection method for solving the % root-finding problem f(x) = 0 for a continuous function f on the % closed interval [a,b] % % Need f(a) and f(b) to have different signs % % Input: f = @(x) function % a = left endpoint of x interval % b = right endpoint of x interval % tol = tolerance for stopping criterion % maxiter = maximum number of iterations % verbose = prints extra information if equal to 1 % % Output: p = approximated root of f on [a,b] % iter = total number of iterations performed % relerr = resulting relative error % p_all = approximated root of f on [a,b] at each step % iter_all = total number of iterations performed at each step % relerr_all = resulting relative error at each step iter=0; FA=f(a); %Initialize output arrays p_all=[]; iter_all=[]; relerr_all=[]; if verbose fprintf("Bisection Method Results:\n"); end while iter=0) || (FP>-1e-15 && FP<=0) || relerr 0 a=p; FA=FP; else b=p; end end if verbose if iter>maxiter fprintf('Method failed after %.0f iterations.\n', maxiter); end end end %% CentralDiff function [df_x0] = CentralDiff(f,h) % Central-Difference Formula for approximating first-derivative % at x0. % % Implemented per Equation 4.1 of Burden, Faires, Burden. % % Input: f = Vector [f(x0+h) f(x0-h)] % h = Distance between x-nodes. % % Output: df_x0 = Approximation of f'(x0) % Apply forward difference approximation df_x0 = 1/(2*h)*(f(1)-f(2)); end %% EstimateFPIter function [N] = EstimateFPIter(g,k,p0,tol,verbose) % Function for estimating the number of iterations required to achieve % the desired tolerance using Fixed Point Iteration. Based on corollary % to Fixed-Point Theorem. % % Input: f = @(x) function % k = constant from Fixed-Point Theorem (0 < k < 1) % p0 = any number on the interval [a,b] % tol = tolerance for stopping criterion % % Output: N = estimated number of iterations to achieve tolerance % Calculate p1 p1=g(p0); %Calculate left hand side of inequality and move terms over left_side=tol*(1-k)/abs(p1-p0); %Use log property: logB(X) = logA(X) / logA(B) N=log(left_side) / log(k); if verbose fprintf('Number of iterations required to converge '); fprintf('to fixed point within tolerance of %.8f is <=%.4f.\n',tol,N); end end %% EulersMethod function [tt,w]=EulersMethod(f,a,b,h,alpha,method) % Uses Euler's Method to approximate the solution of a well-posed IVP % at equally spaced numbers in the interval [a,b]. % % Implemented per p.267 of Burden, Faires, Burden. % % Input: f = anonymous function y'=f(t,y) % a = start point of interval % b = end point of interval % h = step-size between mesh points % alpha = initial condition y(a)=alpha % method = difference method to use ('forward', 'midpoint') % % Output: tt = mesh points % w = approximation to y at each mesh point % Create time-vector (mesh points). tt=a:h:b; % Initialize output vector w=zeros(length(tt),1); w(1)=alpha; switch method case 'forward' % For each meshpoint, use forward approx. for ii=2:length(tt) % Calculate the y-approximation w(ii)=w(ii-1)+h.*f(tt(ii-1),w(ii-1)); end case 'backward' % For each meshpoint, use backward approx. tol=10^-6; maxiter=100; verbose=0; for ii=2:length(tt) % Calculate the y-approximation
g=@(u) w(ii-1)+h.*f(tt(ii),u)-u; % Use Bisection Method since backward-euler is implicit [w(ii),~,~,~,~,~] = Bisection(g,w(ii-1)-1,w(ii-1)+1,... tol,maxiter,verbose); end case 'midpoint' % Get second initial point using forward approx. w(2)=w(1)+h.*f(tt(1),w(1)); % For each meshpoint afterwords, use midpoint method for ii=3:length(tt) % Calculate the y-approximation
w(ii)=w(ii-2)+2.*h.*f(tt(ii-1),w(ii-1)); end otherwise fprintf("Error, unknown method. "); fprintf("Use 'forward' or 'midpoint'\n"); end end %% FivePointMidpoint function [df_x0] = FivePointMidpoint(f,h) % Five-Point Endpoint Formula for approximating first-derivative % at x0. % % Implemented per Equation 4.6 of Burden, Faires, Burden. % % Input: f = Vector [f(x0-2*h) f(x0-h) f(x0+h) f(x0+2*h)] % h = Distance between x-nodes. % % Output: df_x0 = Approximation of f'(x0) df_x0=1/(12*h).*(f(1)-8*f(2)+8*f(3)-f(4)); end %% FixedPointIter function [p,iter,relerr] = FixedPointIter(g,p0,tol,maxiter) % Function implementing fixed-point iteration for g(x) on [a,b] to % find a solution to p=g(p) given an initial approximation p0. % Assumes g(x) has met existence and uniqueness criteria. % % Algorithm described in Numerical Analysis (Burden, Faires, Burden). % % Input: f = @(x) function % p0 = any number on the interval [a,b] % tol = tolerance for stopping criterion % maxiter = maximum number of iterations % % Output: p = approximated unique fixed point of g on [a,b] % iter = total number of iterations performed % relerr = resulting relative error %Initialize values iter=1; relerr=inf; %Perform fixed point iteration while (iter
=maxiter fprintf('Method failed after %d iterations.\n',iter); end end %% ForwardDiff function [df_x0] = ForwardDiff(f,h) % Forward-Difference Formula for approximating first-derivative % at x0. % % Implemented per Equation 4.1 of Burden, Faires, Burden. % % Input: f = Vector [f(x0+h) f(x0)] % h = Distance between x-nodes. % % Output: df_x0 = Approximation of f'(x0) % Apply forward difference approximation df_x0 = 1/h*(f(1)-f(2)); end %% GaussSeidelMethod function [XO,iter,norm, ... XO_all,iter_all,norm_all] = GaussSeidelMethod(A,b,x0,tol,maxiter,verbose) % Function implementing Gauss-Seidel method for iteratively solving a % linear system of equations. % % Implemented per p.461 of Burden, Faires, Burden. % % Input: A = A-matrix (from form Ax=b) % b = b-matrix (from form Ax=b) % x0 = Initial approximation of x (from form Ax=b) % tol = tolerance for stopping criterion % maxiter = maximum number of iterations % verbose = prints extra information if equal to 1 % % Output: XO = approximated solution for x (from form Ax=b) % iter = total number of iterations performed % norm = resulting infinity norm % XO_all = approximated solution for x at each step % iter_all = total number of iterations performed at each step % norm_all = resulting infinity norm at each step % Identify number of equations/unknowns n=length(b); XO=x0; x=zeros(n,1); iter=0; %Initialize output arrays XO_all=[]; iter_all=[]; norm_all=[]; if verbose fprintf("Gauss-Seidel Method Results:\n"); end % Begin iteration while (iter maxiter fprintf('Method failed after %.0f iterations.\n', maxiter); end end end % GramSchmidt function [phi_k,Bk,Ck]=GramSchmidt(w,a,b,k,prev_phi,pprev_phi) % Uses Gram-Schmidt process to construct orthogonal polynomials % on the interval [a,b] % % Input: w = weight function % a = start point of interval % b = end point of interval % k = current k-value % prev_phi = k-1 polynomial function % pprev_phi = k-2 polynomial function % % Output: phi_k = current (k) polynomial function (anonymous) % Bk = B coefficient % Ck = C coefficient % Calculate Bk B_num=@(x) x.*w.*(prev_phi(x).^2); B_den=@(x) w.*(prev_phi(x).^2)+0.*x; Bk=double(integral(B_num,a,b)./integral(B_den,a,b)); if k==1 Ck=0; % Define polynomial function for k==1 phi_k=@(x) x-Bk; else % Calculate Ck for k>2 C_num=@(x) x.*w.*prev_phi(x).*pprev_phi(x); C_den=@(x) w.*(pprev_phi(x).^2)+0.*x; Ck=double(integral(C_num,a,b)./integral(C_den,a,b)); % Define polynomial function for k>2 phi_k=@(x) (x-Bk).*prev_phi(x)-Ck.*pprev_phi(x); end end %% Jacobian function J=Jacobian(F,X) % Function for symbollically calculating the Jacobian for a system of % equations and unknown variables. % % Input: F = Column cell array of symbolic functions % X = Column vector of symbolic unknown variables % % Output: J = Symbolic Jacobian matrix (n x n) % Get number of input functions n=length(F); % Make sure number of functions and unknowns match if n == length(X) % Initialize Jacobian as a symbolic matrix J=sym('j',[n n]); % For each row for ii=1:n % For each column for jj=1:n % Calculate df(ii)/dx(jj) J(ii,jj)=diff(F(ii),X(jj)); end end else disp('Number of functions and unknowns do not match.'); end end %% JacobisMethod function [XO,iter,norm, ... XO_all,iter_all,norm_all] = JacobisMethod(A,b,x0,tol,maxiter,verbose) % Function implementing Jacobi's method for iteratively solving a % linear system of equations. % % Implemented per p.459 of Burden, Faires, Burden. % % Input: A = A-matrix (from form Ax=b) % b = b-matrix (from form Ax=b) % x0 = Initial approximation of x (from form Ax=b) % tol = tolerance for stopping criterion % maxiter = maximum number of iterations % verbose = prints extra information if equal to 1 % % Output: XO = approximated solution for x (from form Ax=b) % iter = total number of iterations performed % norm = resulting infinity norm % XO_all = approximated solution for x at each step % iter_all = total number of iterations performed at each step % norm_all = resulting infinity norm at each step % Identify number of equations/unknowns n=length(b); XO=x0; x=zeros(n,1); iter=0; %Initialize output arrays XO_all=[]; iter_all=[]; norm_all=[]; if verbose fprintf("Jacobi's Method Results:\n"); end % Begin iteration while (iter maxiter fprintf('Method failed after %.0f iterations.\n', maxiter); end end end %% Lagrange function poly = Lagrange(xpts,ypts,xeval) % Function to generate Lagrange interpolating polynomial at % values xeval that passes through points (xpts,ypts) % % Input: xpts = x points % ypts = y points % xeval = evaluate polynomial at these x values % % Output: poly = Lagrange interpolating polynomial (order n) N = length(xpts); %n+1 L = ones(N,length(xeval)); poly = 0; % Generate Lagrange functions L_k(x) for each point xeval for k = 1:N for i = 1:N if (i ~= k) L(k,:) = L(k,:).*(xeval-xpts(i))./(xpts(k)-xpts(i)); end end % Generate Lagrange interpolating polynomial poly = poly + ypts(k)*L(k,:); end end %% LLS function [theta]=LLS(x,y,deg) % Function for constructing the linear least squares (LLS) polynomial % coefficients % % Input: x = Input datapoints/nodes, xn, column vector % y = Input datapoints, yn=f(xn), column vector % deg = Degree of polynomial % % Output: theta = linear least squares polynomial coefficients % Construct design matrix X = zeros(length(x),deg); for ii=0:deg X(:,ii+1)=x.^ii; end % Solve for coefficients theta = (X'*X)\(X'*y); % Reverse order of theta for use with polyval (descending powers) theta=flip(theta); end %% LogB function [logBX] = LogB(X,B) % Function for calculating log(x) with any base % using log property: logB(X) = logA(X) / logA(B) % % Input: X = value to take log of % B = Base % % Output: logBX = log(X) with base 'B' logBX = log(X)/log(B); end %% NaturalCubicSpline function [a,b,c,d,S] = NaturalCubicSpline(x,y) % Function for constructing the cubic spline interpolant S for the % function f, defined at x0 maxiter fprintf('Method failed after %.0f iterations', maxiter); end end %% NewtonsMethodForSystems function [XO,iter,norm, ... XO_all,iter_all,norm_all] = NewtonsMethodForSystems(F,X,x0,tol,maxiter,verbose) % Function implementing Newton's Method for iteratively solving a % nonlinear system of equations F(x)=0. % % Implemented per p.653 of Burden, Faires, Burden. % % Input: F = Column cell array of nonlinear mappings % X = Column vector of symbolic variables % x0 = Initial approximation of x % tol = tolerance for stopping criterion % maxiter = maximum number of iterations % verbose = prints extra information if equal to 1 % % Output: XO = approximated solution for x (from form Ax=b) % iter = total number of iterations performed % norm = resulting infinity norm % XO_all = approximated solution for x at each step % iter_all = total number of iterations performed at each step % norm_all = resulting infinity norm at each step % Initalize variables XO=x0; iter=0; % Initialize output arrays XO_all=[]; iter_all=[]; norm_all=[]; if verbose fprintf("Newton's Method for Systems Results:\n"); end % Begin iteration while (iter maxiter fprintf('Method failed after %.0f iterations.\n', maxiter); end end end %% QuasiNewton function [p,iter,relerr, ... p_all,iter_all, ... relerr_all]=QuasiNewton(f,a,b,tol,max_global_steps, ... max_local_steps,verbose) % Function implementing Quasi-Newton scheme for solving the root-finding % problem f(x) = 0 for a continuous function f on the closed interval % [a,b]. Uses Bisection as global method to approximate p0, then % Newton's Method for final convergence. % % Input: f = @(x) function % a = left endpoint of x interval % b = right endpoint of x interval % tol = tolerance for stopping criterion % max_global_steps = maximum number of Bisection iterations % max_local_steps = maximum number of Newton's Method iterations % verbose = prints extra information if equal to 1 % % Output: p = approximated root of f on [a,b] % iter = total number of iterations performed % relerr = resulting relative error % p_all = approximated root of f on [a,b] at each step % iter_all = total number of iterations performed at each step % relerr_all = resulting relative error at each step % iter_type = identify whether global or local step % Perform global bisection steps to get a good p0 approximation [p0,iter_global,relerr_global, ... p0_all,iter_global_all, ... relerr_global_all] = Bisection(f,a,b,tol,max_global_steps,0); % Store global results in final output arrays p_all=p0_all; iter_all=iter_global_all; relerr_all=relerr_global_all; iter_type(1:length(p0_all))="global"; % Perform local Newton's Method steps until convergence [p,iter,relerr, ... p_local_all,iter_local_all, ... relerr_local_all] = NewtonsMethod(f,p0,tol,max_local_steps,0); % Add local results to final output arrays iter_local_all=iter_local_all+iter_global; p_all=[p_all,p_local_all]; iter_all=[iter_all,iter_local_all]; relerr_all=[relerr_all,relerr_local_all]; iter_type(length(p0_all)+1:length(p_all))="local"; if verbose fprintf("Quasi-Newton Results:\n"); for ii=1:length(p_all) fprintf('iter: %.f, ',iter_all(ii)); fprintf('p%.f: %.8f, ',iter_all(ii),p_all(ii)); fprintf('relerr: %.4f, ', relerr_all(ii)); fprintf('type: %s\n', iter_type(ii)); end end end %% QuasiSecant function [p,iter,relerr, ... p_all,iter_all, ... relerr_all]=QuasiSecant(f,a,b,tol,max_global_steps, ... max_local_steps,verbose) % Function implementing Quasi-Secant scheme for solving the root-finding % problem f(x) = 0 for a continuous function f on the closed interval % [a,b]. Uses Bisection as global method to approximate p0 and p1, then % Secant Method for final convergence. % % Input: f = @(x) function % a = left endpoint of x interval % b = right endpoint of x interval % tol = tolerance for stopping criterion % max_global_steps = maximum number of Bisection iterations % max_local_steps = maximum number of Secant Method iterations % verbose = prints extra information if equal to 1 % % Output: p = approximated root of f on [a,b] % iter = total number of iterations performed % relerr = resulting relative error % p_all = approximated root of f on [a,b] at each step % iter_all = total number of iterations performed at each step % relerr_all = resulting relative error at each step % iter_type = identify whether global or local step % Perform global bisection steps to get a good p0 and p1 approximation [p_global,iter_global,relerr_global, ... p_global_all,iter_global_all, ... relerr_global_all] = Bisection(f,a,b,tol,max_global_steps,0); % Index results for p0 and p1 p0=p_global_all(length(p_global_all)-1); p1=p_global_all(length(p_global_all)); % Store global results in final output arrays p_all=p_global_all; iter_all=iter_global_all; relerr_all=relerr_global_all; iter_type(1:length(p_global_all))="global"; % Perform local Secant Method steps until convergence [p,iter,relerr, ... p_local_all,iter_local_all, ... relerr_local_all] = SecantMethod(f,p0,p1,tol, ... max_local_steps,0); % Add local results to final output arrays iter_local_all=iter_local_all+iter_global; p_all=[p_all,p_local_all]; iter_all=[iter_all,iter_local_all]; relerr_all=[relerr_all,relerr_local_all]; iter_type(length(p_global_all)+1:length(p_all))="local"; if verbose fprintf("Quasi-Secant Results:\n"); for ii=1:length(p_all) fprintf('iter: %.f, ',iter_all(ii)); fprintf('p%.f: %.8f, ',iter_all(ii),p_all(ii)); fprintf('relerr: %.4f, ', relerr_all(ii)); fprintf('type: %s\n', iter_type(ii)); end end end %% SecantMethod function [p,iter,relerr, ... p_all,iter_all,relerr_all] = SecantMethod(f,p0,p1,tol, ... maxiter,verbose) % Function implementing the Secant method for solving the root-finding % problem f(x) = 0 given initial approximations p0 and p1. % % Implemented per p.71 of Numerical Analysis by Burden, Faires, Burden. % % Input: f = @(x) function % p0 = Initial approximation to p in [a,b] % p1 = Initial approximation of p1 % tol = tolerance for stopping criterion % maxiter = maximum number of iterations % verbose = prints extra information if equal to 1 % % Output: p = approximated root of f % iter = total number of iterations performed % relerr = resulting relative error % p_all = approximated root of f on [a,b] at each step % iter_all = total number of iterations performed at each step % relerr_all = resulting relative error at each step % Perform Secant Method iter=1; q0=f(p0); q1=f(p1); if verbose fprintf("Secant Method Results:\n"); fprintf('iter: 0, p0: %.16f\n', p0); fprintf('iter: 1, p1: %.16f\n', p1); end %Initialize output arrays p_all=[p1]; iter_all=[1]; relerr_all=[abs(p1-p0)]; while iter maxiter fprintf('Method failed after %.0f iterations', maxiter); end end %% RK2 function y = RK2(f,t,y0) % Function to compute RK2 approximation of solution to IVP y'=f(t,y) with % initial value y0 for t over the interval [a,b] % % Input: f = RHS function of DE % t = mesh points of t values over the interval [a,b] % y0 = initial value % % Output: y = RK2 approximation of solution to IVP h = t(end)-t(end-1); % step size N = length(t); % number of points in mesh / approximation y = NaN(N,1); y(1) = y0; for i = 1:N-1 t_half = t(i)+h/2; % time points halway between mesh points y(i+1) = y(i) + h*f([t_half,y(i)+h/2*f([t(i),y(i)])]); % RK2 formula end end %% TayloyPoly function [Pn]=TaylorPoly(f,x0,N) % Function for symbolically calculating the first N-terms of a % Taylor Polynomial. % % Input: f = anonymous @(x) function to approximate % x0 = centerpoint of Taylor series % N = number of terms to calculate % % Output: Pn = symbolic first N-terms of Taylor Polynomial % w.r.t. x. syms x; Pn=0; for k=0:N %Get kth derivative of f fk=diff(f,x,k); %Calculate current term of polynomial Pn_k=subs(fk,x0)/factorial(k)*((x-x0)^k); %Add current term to total Pn=Pn+Pn_k; end end %% TaylorPolyNTerm function [PN]=TaylorPolyNTerm(f,x0,N) % Function for symbolically calculating the Nth term only of a % Taylor Polynomial. % % Input: f = anonymous @(x) function to approximate % x0 = centerpoint of Taylor series % N = number of term to calculate % % Output: PN = symbolic Nth term of Taylor Polynomial w.r.t. x. syms x; %Get Nth derivative of f fN=diff(f,x,N); %Calculate current term of polynomial PN=subs(fN,x0)/factorial(N)*((x-x0)^N); end %% ThreePointEndpoint function [df_x0] = ThreePointEndpoint(f,h) % Three-Point Endpoint Formula for approximating first-derivative % at x0. % % Implemented per Equation 4.4 of Burden, Faires, Burden. % % Input: f = Vector [f(x0) f(x0+h) f(x0+2*h)] % h = Distance between x-nodes. % % Output: df_x0 = Approximation of f'(x0) df_x0=1/(2*h).*(-3*f(1)+4*f(2)-f(3)); end %% ThreePointMidpoint function [df_x0] = ThreePointMidpoint(f,h) % Three-Point Midpoint Formula for approximating first-derivative % at x0. % % Implemented per Equation 4.5 of Burden, Faires, Burden. % % Input: f = Vector [f(x0+h) f(x0-h)] % h = Distance between x-nodes. % % Output: df_x0 = Approximation of f'(x0) df_x0=1/(2*h).*(f(1)-f(2)); end %% TruncationError function [Rn]=TruncationError(f,x0,n,ksi_x) % Function for symbolically calculating the truncation error % associated with Taylor Polynomial approximation. % % Input: f = anonymous @(x) function to approximate % x0 = centerpoint of Taylor series % n = order of Taylor Polynomial % ksi_x = unknown function between x0 and x % % Output: Rn = symbolic Truncation Error w.r.t. x. syms x; %Get (n+1)th derivative of f f_n1=diff(f,x,n+1); %Calculate truncation error Rn=subs(f_n1,ksi_x)/factorial(n+1)*((x-x0)^(n+1)); end %% TruncationErrorLagrange function [Rn]=TruncationErrorLagrange(f,xn,ksi_x) % Function for symbolically calculating the truncation error % associated with Lagrange Interpolating Polynomial approximation. % % Input: f = anonymous @(x) function to approximate % xn = vector of x-points used for interpolation % ksi_x = value of unknown function (between x0 and xn) % % Output: Rn = symbolic Truncation Error w.r.t. x. syms x; %Get (n+1)th derivative of f n=length(xn); f_n1=diff(f,x,n+1); %Calculate truncation error Rn=subs(f_n1,ksi_x)/factorial(n+1); for ii=1:n Rn=Rn*(x-xn(ii)); end end end end
References
- Numerical Analysis, 10th Edition by Burden, Faires, and Burden