Coding Project 1 (Numerical Analysis class , Matlab) Given A function y=f(x), X=
ID: 670312 • Letter: C
Question
Coding Project 1 (Numerical Analysis class , Matlab)
Given A function y=f(x), X=(x0,x1,....,xn)^(T) with yi=f(xi) and x0<x1<....<xn, and one value a, write a Matlab function [y]=//P1(X,n,a) to do the following.
1. Calculate the coefficients of the interpolation polynomial Pn(x) for the points (x0,y0),...,(xn,yn) using Newton divided-differences as in Algorithm 3.2.
2. For the given value a, calculate y= Pn(a). y is the output of your function.
3. Draw the graphs of y= f(x) and y=Pn(x) on the interval [x0,xn] togethert on one coordinate system.
Explanation / Answer
1.
function [Q,fcount] = quadtx(F,a,b,tol,varargin)
%QUADTX Evaluate definite integral numerically.
% Q = QUADTX(F,A,B) approximates the integral of F(x)
% from A to B to within a tolerance of 1.e-6.
%
% Q = QUADTX(F,A,B,tol) uses tol instead of 1.e-6.
%
% The first argument, F, is a function handle or
% an anonymous function that defines F(x).
%
% Arguments beyond the first four,
% Q = QUADTX(F,a,b,tol,p1,p2,...), are passed on to the
% integrand, F(x,p1,p2,..).
%
% [Q,fcount] = QUADTX(F,...) also counts the number of
% evaluations of F(x).
%
% See also QUAD, QUADL, DBLQUAD, QUADGUI.
% Default tolerance
if nargin < 4 | isempty(tol)
tol = 1.e-6;
end
% Initialization
c = (a + b)/2;
fa = F(a,varargin{:});
fc = F(c,varargin{:});
fb = F(b,varargin{:});
% Recursive call
[Q,k] = quadtxstep(F, a, b, tol, fa, fc, fb, varargin{:});
fcount = k + 3;
2.
function a = InterpV(x,y)
% a = InterpV(x,y)
% This computes the Vandermonde polynomial interpolant where
% x is a column n-vector with distinct components and y is a
% column n-vector.
%
% a is a column n-vector with the property that if
%
% p(x) = a(1) + a(2)x + ... a(n)x^(n-1)
% then
% p(x(i)) = y(i), i=1:n
n = length(x);
V = ones(n,n);
for j=2:n
% Set up column j.
V(:,j) = x.*V(:,j-1);
end
a = Vy;
3.
function [c, err, yc] = regula (f, a, b, delta, epsilon, maxl)
%Input - f is the function input as a string ’f’
% - a and b are the left and right end points
% - delta is the tolerance for the zero
% - epsilon is the tolerance for the value of f at the zero
% - max1 is the maximum number of iterations
%Output - c is the zero
% - yc = f(c)
% - err is the error estimate for c
ya=feva1(f,a);
yb=feva1(f,a);
if ya*yb>0
disp (’Note: f(a)*f(b)>0’),
break,
end
for k=1:max1
dx=yb*(b-a)/(yb-ba);
c=b-dx;
ac=c-a;
yc=feval(f,c);
if yc==0, break;
elseif yb*yc>0
b=c;
yb=yc;
end
dx=min (abs (dx), ac);
if abs (dx)<delta, break, end
if abs (yc)<epsilon, break, end
end
c;
err=abs (b-a)/2;
yc=feval(f,c);