Mathematics Homework Help

https://www.youtube.com/watch?time_continue=1058&v=oENK6VIgbqU&feature=emb_title video is related to the question 7

Textbook Codes MYSPLINE and MYSPLINE_VAL

function y = myspline_val(pp, x, deriv)

% Evaluate the cubic spline calculated in myspline.m or any

% derivatives at x

%%% input:

% pp = structure containing coefficients and endpoints of subintervals

% which can be the output of interp1 or myspline or csa

% x = points at which cubic spline is evaluated

% deriv = (optional)

% = (omitted) derive = 0

% >= 0 the value or the derivative which is desired

%%% output:

% y = pp(x) or pp'(x) or pp”(x) or …

if nargin == 2

deriv = 0;

end

lenx1 = size(x, 1);

x = x(:);

[~, ~, ixdp] = histcounts(x, [-inf, pp.breaks(2:end-1), +inf]);

C = fliplr(pp.coefs);

if deriv == 1

C = C*diag([1 2 3], -1);

elseif deriv == 2

C = C*diag([2 6], -2);

elseif deriv == 3

C = C*diag([6], -3);

elseif deriv > 3

C = zeros(size(C));

end

xx = x – pp.breaks(ixdp)’;

y = C(ixdp,1) + xx.*(C(ixdp,2) + xx.*(C(ixdp,3) + xx.*C(ixdp,4)));

if lenx1 == 1

y = y.’; % uses the transpose, not the conjugate transpose

end

end

function pp = myspline(xdp, ydp, bc)

% MYSPLINE Generate a cubic spline given the data points

% {(xdp_i,ydp_i) | i = 1, 2, …, n}

% The boundary conditions are contained in the structure bc .

% They can be different on the left end, i.e., xdp(1) , and the

% right end, i.e., xdp(end)

%%% input:

% xdp, ydp = the data points

% bc = (optional, default: use not-a-knot for both left and right boundaries)

% = structure containing the boundary conditions

% = [] – same as omitted

% .l_method = ‘1’ – first derivative input, value is in bc.l_val

% ‘2’ – second derivative input, value is in bc.l_val

% ‘3’ – third derivative input, value is in bc.l_val

% ‘notaknot’ – not-a-knot boundary condition

% ‘periodic’ – periodic boundary conditions

% .r_method = ‘1’ – first derivative input, value is in bc.l_val

% ‘2’ – second derivative input, value is in bc.l_val

% ‘3’ – third derivative input, value is in bc.l_val

% ‘notaknot’ – not-a-knot boundary condition

% NOTE: if the left boundary condition is periodic, the

% right must be also, so it is not needed (but you

% are free to set it to ‘periodic’ if you wish)

% HOWEVER, if the right one has been assigned

% anything other than ‘periodic’, a fatal error message is

% returned

%%% output:

% pp = structure containing spline data which can be input directly to myspline_val

% .breaks = a row vector of breakpoints

% .coefs = (n-1)x4 matrix whose i-th row contains coefficients of p_i(x)

% in descending order when written in shifted power form

%

% C_{i,1} (x-x_i)^3 + C_{i,2} (x-x_i)^2 + C_{i,3} (x-x_i) + C_{i,4}

%

% Originally from “Learning MATLAB, Problem Solving, and Numerical Analysis

% Through Examples” by Overman

% Modified by Tae Eun Kim, July 2020.

%% If there are only two inputs or if bc is []

if nargin == 2 || isempty(bc)

bc.l_method = ‘notaknot’;

bc.r_method = ‘notaknot’;

end

%% Defining key expressions

n = length(xdp);

% dx_i = x_{i+1} – x_i for i = 1, 2, …, n-1

% dy_i = y_{i+1} – y_i for i = 1, 2, …, n-1

% dydx_i = dy_i/dx_i for i = 1, 2, …, n-1

% nx_i = dx_{i-1} + dx_i for i = 2, …, n-1

% NOTE: nx_1 is needed for periodic bc

xdp = xdp(:);

ydp = ydp(:);

dx = diff(xdp);

dy = diff(ydp);

dydx = dy./dx;

nx = zeros(n-1, 1);

nx(2:n-1) = dx(1:n-2) + dx(2:n-1);

%% Setting up X * sigma = r

% X = nxn matrix

% initially fill in everything but the first and last rows

X = 2*diag([0; nx(2:n-1); 0]) + …

diag([0; dx(1:n-2)], 1) + …

diag([dx(2:n-1); 0], -1);

% r = nx1 vector of the right-hand side

% initially fill in everything but the first and last elements

r = zeros(n, 1);

r(2:n-1) = 3*( dydx(1:n-2).*dx(2:n-1) + …

dydx(2:n-1).*dx(1:n-2) );

%% Boundary Conditions

X(1,1:3) = [dx(2)^2, dx(2)^2 – dx(1)^2, -dx(1)^2];

r(1) = 2*( dydx(1)*dx(2)^2 – dydx(2)*dx(1)^2 );

X(n,n-2:n) = [dx(n-1)^2, dx(n-2)^2 – dx(n-1)^2, dx(n-2)^2];

r(n) = 2*( dydx(n-1)*dx(n-2)^2 – dydx(n-2)*dx(n-1)^2 );

%% Solving X * sigma = r

sigma = sparse(X)r;

%% Packaging Outputs

pp.breaks = xdp’;

pp.coefs = [( sigma(1:n-1) + sigma(2:n) – 2*dydx )./( dx.^2 ), …

( 3*dydx – 2*sigma(1:n-1) – sigma(2:n) )./( dx ), …

sigma(1:n-1), …

ydp(1:n-1)];

end

Usage

Clamped

xdp = 0:5;

f = @(x) x.^3;

ydp = f(xdp);

bc.l_method = ‘1’;

bc.r_method = ‘1’;

bc.l_val = 0;

bc.r_val = 75;

x = linspace(xdp(1), xdp(end), 10000);

pp = myspline(xdp, ydp, bc);

y = myspline_val(pp, x);

plot(x, y, xdp, ydp, ‘ro’, x, f(x), ‘:’);

norm(y-f(x), ‘Inf’)