• 对不起，我没有MatLab，但这是一个类似的数学软件。我觉得你可以把它翻译出来*

与answer of Savithru中所示的倾斜圆柱体不同，我们可以使用遵循给定曲线的版本。
从本质上讲，你需要的是一组圆，它们都与你试图沿着的曲线f(x)垂直。首先，让我们定义一个圆：

圈子：

假设在垂直于单位向量w=(w1,w2,w3)的平面上有一个半径为R的圆，它穿过点(X0,Y0,Z0)。该平面由两个垂直于w的单位向量u=(u1,u2,u3)和v=(v1,v2,v3)定义。圆的参数方程为：

x = X0 + R cos(theta) u1 + R sin(theta) v1
y = Y0 + R cos(theta) u2 + R sin(theta) v2
z = Z0 + R cos(theta) u3 + R sin(theta) v3

其中，theta在0到2π的间隔内变化。
现在我们定义了圆，让我们定义管子。

管子：

要定义曲面，我们需要两个参数，第一个参数是圆的theta。如果我们想跟踪f(x)，第二个将是x。我们所需要做的就是定义正交向量u和v。我们知道它们与单位向量w是正交的，这只不过是用一阶导数得到的f(x)的切线。因此，您可以定义：

w = Normalize[{1,f'(x),0}]
u = Normalize[Cross[w,{0,0,1}]]
v = Cross[w,u]

因此，您的参数方程变为：

x = x    + R cos(theta) u1(x) + R sin(theta) v1(x)
y = f(x) + R cos(theta) u2(x) + R sin(theta) v2(x)
z = 0    + R cos(theta) u3(x) + R sin(theta) v3(x)

在MATHEMICAL A中，这是这样的：

R=1
f[x_] := Sin[x]
w[x_] := Normalize[{1, f'[x], 0}]
u[x_] := Normalize[Cross[w[x], {0, 0, 1}]]
v[x_] := Cross[w[x], u[x]]
ParametricPlot3D[{x, f[x], 0} + R Cos[t] u[x] + R Sin[t] v[x], {x, 0, 2 Pi}, {t, 0, 2 Pi}]

**注意：**您可以使用从3D空间f(r)中的曲线的切线、法线和双法线向量构建的Frenet-Serret帧轻松扩展这一点

首先，应指定圆柱体的轴的方向。现在，我假设它指向z方向，它将只在x方向振荡(即轴的方程是x = sin(z)和y=0)。
如果圆柱体的轴随z变化，则圆柱体表面的x,y坐标也应该是z的函数。您可以这样做：首先计算直线圆柱体的x,y点，然后添加一个依赖于本地z值的“Shift”。
代码如下：

% Parameters
r=5; l=5; nTheta=100, nL = 20;

theta = linspace(0,2*pi,nTheta+1);
x = r * cos(theta);
y = r * sin(theta);

z = linspace(0,l,nL)';
xshift = repmat( sin(z), 1, nTheta+1); %this is a function of z

X = repmat(x,nL,1) + xshift;
Y = repmat(y,nL,1);
Z = repmat(z, 1, nTheta+1);

% Plot
surf(X,Y,Z)

如果圆柱体的轴在x方向和y方向上都振荡(或曲线)，则还需要yshift。

我想我应该在这里加上我的答案，因为我正试图做类似的事情，而@kvantour的答案真的很有帮助。在我的例子中，我使用了Python，并将圆柱体的两端逐渐变细。

import numpy as np

length = 1  # Length of the cylinder
N = 128  # Number of vertices along midline
N_theta = 32  # Number of rotation angles
r_max = 0.02  # Maximum radius
N_taper = int(N / 4)  # Number of points to taper at the ends
shape_k1 = 1.5  # Changes shape at midpoint
shape_k2 = 1  # Changes shape at endpoint
u = np.linspace(0, length, N)
theta = np.linspace(0, 2 * np.pi, N_theta)

# Midline as straight line
X = np.stack([u, np.zeros_like(u), np.zeros_like(u)], axis=-1)

# Midline as a spiral
X = np.zeros((N, 3))
X[:, 0] = np.sin(2 * np.pi * u) / 10
X[:, 1] = np.cos(2 * np.pi * u) / 10
X[:, 2] = np.linspace(1 / np.sqrt(3), 0, N)

# Define the radius as a function of midline
s = np.linspace(0, np.pi / 2, N_taper)
x = N_taper * np.cos(s)**shape_k1
y = r_max * np.sin(s)**shape_k2
t_idxs = np.linspace(0, N_taper, N_taper)

# Resample to get equally spaced points
f = interp1d(x, y, kind='cubic', fill_value='extrapolate')
slopes = f(t_idxs)
r = np.concatenate([
    slopes[::-1],
    np.ones(N - 2 * N_taper) * r_max,
    slopes
])

# Calculate a TNB frame from the midline
T = normalise(np.gradient(X, axis=0))
N = normalise(np.cross(T, np.ones_like(T) * [0, 0, 1]))
B = np.cross(T, N)

# Generate surface
x = X[:, 0] + r * (np.outer(np.cos(theta), N[:, 0]) + np.outer(np.sin(theta), B[:, 0]))
y = X[:, 1] + r * (np.outer(np.cos(theta), N[:, 1]) + np.outer(np.sin(theta), B[:, 1]))
z = X[:, 2] + r * (np.outer(np.cos(theta), N[:, 2]) + np.outer(np.sin(theta), B[:, 2]))

# (Also equivalent to):
surface = X + r[None, :, None] \
              * (np.einsum('i,jk->ijk', np.cos(theta), N) + np.einsum('i,jk->ijk', np.sin(theta), B))
x, y, z = surface[..., 0], surface[..., 1], surface[..., 2]

# Plot using matplotlib (doesn't look great without tweaking axis limits)
import matplotlib.pyplot as plt
ax = plt.axes(projection='3d')
ax.plot_surface(x, y, z, cmap='viridis')
plt.show()

# Plot using mayavi - looks much nicer out the box
from mayavi import mlab
mlab.mesh(x, y, z)
mlab.show()

Matplotlib：

Mayavi：