Suppose $f'(x_1),\ f'(x_2),\ f'(x_3)$ are given, how to give a polynomial interpolation $p(x)$ such that $p'(x_1)= f'(x_1),\ p'(x_2)=f'(x_2),\ p'(x_3)=f'(x_3)$? And how to give an error analysis?
The problem boils down to the solution of a linear system a equations. You get one equation for the function value at $f(x_1)$, and three more equations from your derivative equations. Then, a polynomial of degree three is suitable.
A more brute-force, not so clever way is to use computer algebra such as SymPy for a symbolical solution:
from sympy import * N = 3 # degree of polynomial - for N + 1 'pieces of information', use a polynomial of degree N known_x = symbols(["x_%i"%i for i in range(N)]) # create symbols for supporting points known_derivs = symbols(["d_%i"%i for i in range(N)]) # create symbols for supporting points y1 = Symbol("y_1") # the single zeroth derivative we know # create symbols for the independent variable and the coefficients: x = Symbol("x") a = [ Symbol("a_%i"%i) for i in range(N+1) ] # generate polynomial and the derivative p = sum([ coeff*x**i for i, coeff in enumerate(a) ]) p_prime = diff(p, x, 1) # define the equations to solve: eqns = [ p_prime.subs(x, xi) - d for xi, d in zip(known_x, known_derivs) ] # equation for the first derivative eqns.append(p.subs(x, known_x) - y1) # solve for the coefficients and print solution: solution = solve(eqns, a) pretty_sol = [ str(coeff) + " = " + str(value) + "\n" for coeff, value in zip(solution.keys(), solution.values()) ] print("The coefficients are:\n" + "".join(pretty_sol))
This not very elegant, but seems to work even you added another equation for the derivative (so that could use a polynomial of degree 'N = 4' for the five equations).