Is it possible to do an integration over a Bézier curve "area"? I have a shape defined by a Bézier spline that has a width,
and I want to give it an attractive force.
Is such a thing even doable without approximating it?

Edit:
To integrate the area, I think I will first have to express the answer for a line (crossection) of the path, then integrate it through the path.

Edit:
Moving into the teritory of actually answering the question, but I think I'll keep editing the question until I actually have a final answer.
So, now I want to rename the spline, it will new be defined by $$Q(t) = (1-t)^3\ Q_0+(1-t)^2t\ Q_1+(1-t)t^2\ Q_2+t^3\ Q_3$$ rearranging gives $$Q_3\ t^3+Q_2\ (t^2-t^3)+Q_1(t-2t^2+t^3)+Q_0(1-3t+3t^2-t^3)\\=\\(\underbrace{Q_3-Q_2+Q_1-Q_0}_{\mathcal{A}})\ t^3+(\underbrace{Q_2-2Q_1+3Q_0}_{\mathcal{B}})\ t^2+(\underbrace{Q_1-3Q_0}_{\mathcal{C}})\ t+Q_0$$
Finding the speed vector and its square $Q′(t)$:$$Q′(t)=3\mathcal{A}t^2+2\mathcal{B}t+\mathcal{C}$$
$$(Q′(t))^2=9\mathcal{A}^2t^4+4\mathcal{B}^2t^2+\mathcal{C}^2+12\mathcal{A}\mathcal{B}t^3+6\mathcal{A}\mathcal{C}t^2+4\mathcal{B}\mathcal{C}t$$
Now I'm going to find $P_0(t)$ and $P_1(t)$.

Then, I want to rotate the speed vector left and right, putting $P_0$ on the left side and $P_1$ on the right. Rotating a vector $(x,\ y)$ gives $(-y,\ x)$. Rotating right gives $(y,\ -x)$.
We also devide by the length of the speed vector to get a unit vector. Then multiplying with the spline width $W$.
We name the length of the speed vector $$\Gamma(t)=\sqrt{Q′_x(t)^2+Q′_y(t)^2}$$
We get $$P_0(t) = Q(t) + W{(-Q′_y(t),\ Q′_x(t))\over\Gamma(t)}\\[1.5em]P_1(t) = Q(t) + W{(Q′_y(t),\ -Q′_x(t))\over\Gamma(t)}$$
Can I now just do this? $$\iint{1 \over \lVert P(u,t)\rVert ^2}\Gamma(t)\ du\ dt$$
I'm trying to express: the force of attraction by a spline cross section times how far the cross section moves at time $t$.
 A: Take a look at the following illustration featuring a curve similar to the one you have shown and the acceleration vectors that have the nice (ill-known) representation:
$$a=(1-t)(A-2B+C)+t(B-2C+D) \ \ \text{for} \ \ t \in [0,1]$$
Isn't it a tool for expressing what you call attraction ?

Matlab program for the figure:
clear all;close all,hold on;
P=[(4-2*i) (5-i) (5+i) (3+2i)]; % points A,B,C,D
plot(P,'b');
t=0:0.01:1;s=1-t;
m=P(1)*s.^3+P(2)*3*s.^2.*t+P(3)*3*t.^2.*s+P(4)*t.^3;
plot(m,'r');
a=(1-t)*(P(1)-2*P(2)+P(3))+t*(P(2)-2*P(3)+P(4));
quiver(real(m),imag(m),real(a),imag(a),1); % "acceleration field"

A: Reposting an answer that I got on paper for posterity.
Part 1
$$C(t)=P_0B_{0,3}+P_1B_{1,3}+P_2B_{2,3}+P_3B_{3,3}$$ where B is a Bernstein polynomial
$$t\in[0,1]$$
Force on the line:
$$F_{line}=\int_{t=0}^1{1\over(Q(t))^2}\ dt\\[2em]
=\int_{t=0}^1{1\over((1-t)^3\ Q_0+(1-t)^2t\ Q_1+(1-t)t^2\ Q_2+t^3\ Q_3))^2}\ dt$$
Assume: Roots $\alpha$, $\beta$, $\gamma$ of $Q(t)$ exist
$${1\over((t-\alpha)(t-\beta)(t-\gamma))^2}\\[2.5em]=\overbrace{{A\over(t-\alpha)}+{B\over(t-\beta)}+{C\over(t-\gamma)}+{D\over(t-\alpha)^2}+{E\over(t-\beta)^2}+{F\over(t-\gamma)^2}}^{constants}\\[3.5em]
\begin{align}\int_0^1{A\over t-\alpha}\ dt &= A(log(1-\alpha)-log(-\alpha))\\[1.5em]
\int_0^1{D\over (t-\alpha)^2}\ dt &= \left[\begin{array}{}u=t-\alpha\\du=dt\\t=0\implies u=-\alpha\\t=1\implies u=1-\alpha\end{array}\right]\\[1.5em]&=D\int_{-\alpha}^{1-\alpha}{1\over u^2}\ du=D\left[{-1\over u}\right]_{-\alpha}^{1-\alpha}\\[1.5em]&=D\left({-1\over 1-\alpha}+{1\over\alpha}\right)\end{align}$$
So $\int_0^1{1\over(Q(t))^2}\ dt$ turns into six solvable integrals.
As an approximation of your desired force over the area, you could do $F_{area}\approx wF_{line}$. $w$ is width of band.

Part 2
More serious attempt. No approximations.

Say we have two equations $Q_l(t)$ and $Q_r(t)$ that fulfill $${Q_l(t)+Q_r(t)\over2}=Q(t)$$
$$Q_l(t) = Q(t) + \delta\ \vec{v}(t)$$
$$Q_r(t) = Q(t) - \delta\ \vec{v}(t)$$
$\vec{v}(t)$ is always orthogonal to $Q(t)$ and has length 1
$$\vec{v}(t)={\vec{T}′(t)\over\lVert\vec{T}′(t)\rVert},\ \text{where}\ \vec{T}(t)={\vec{Q}′(t)\over\lVert\vec{Q}′(t)\rVert}$$
So $Q_l(t)$ and $Q_r(t)$ exist and should be not too hard to find. We could also construct them by moving $Q$'s control points.
Now introduce a new variable $\omega$ that goes from 0 to 1, and swap out $Q(t)$ in the original  problem with $(1-\omega)\ Q_l(t)+\omega\ Q_r(t)$. We can set up the integral $$\int_{\omega=0}^1\int_{t=0}^1{1\over((1-\omega)\ Q_l(t)+\omega\ Q_r(t))^2}\ dt\ d\omega$$
Solve for the $\omega$ integral first
$$\begin{align}\int_{\omega=0}^1{1\over((Q_r-Q_l)\ \omega+Q_l)^2}\ d\omega&=\int_0^1{1\over((Q_r-Q_l)(\omega+{Q_l\over{Q_r-Q_l}}))^2}\ d\omega\\[2.5em]&={1\over(Q_r-Q_l)^2}\int_0^1{1\over(\omega+{Q_l\over Q_r-Q_l})^2}\ d\omega\\[2.5em]&= {1\over(Q_r-Q_l)^2}\left[{-1\over(1+{Q_l\over Q_r-Q_l}}+{Q_r-Q_l\over Q_l}\right]\\[2.5em]&={-1\over((Q_r-Q_l)^2+Q_l(Q_r-Q_l))}+{1\over Q_l(Q_r-Q_l)}\\[2.5em]&={-1\over Q_r^2-Q_rQ_l}+{1\over Q_lQ_r-Q_l^2}\\[2.5em]&={1\over Q_rQ_l-Q_r^2}+{1\over Q_lQ_r-Q_l^2}\end{align}$$
Two very similar terms. Let's solve the integral for one of them, the other will be analogous. $$\require{cancel}\int_{t=0}^1{1\over Q_lQ_r-Q_l^2}\ dt\\[2.5em]
{1\over Q_lQ_r-Q_l^2}={1\over(Q(t)+\delta\ v(t))\ (Q(t)-\delta\ v(t))-(Q(t)+\delta\ v(t))^2}\\[2.5em]{1\over\cancel{Q^2}-(\delta\ v)^2-\cancel{Q^2}-2Q\delta\ v- (\delta\ v)^2}\\[2.5em]={-1\over2(Q\delta\ v+(\delta\ v)^2)}$$
So you get $$-{1\over2\delta}\int_0^1{1\over Q(t)v(t)+\delta\ v(t)^2}\ dt$$
$v(t)$ should be 1st order polynomial.
Denominator is then 4th order polynomial.
Its roots can be found with some closed form equation.
$$-{1\over2\delta}\int_0^1{1\over(t-\epsilon)(t-\zeta)(t-\eta)(t-\theta)}\ dt$$
The integrand can be split into partial fractions, and we get
$$-{1\over2\delta}\left(\int_0^1\overbrace{{G\over(t-\epsilon)}\ dt+\int_0^1{H\over(t-\zeta)}\ dt+\int_0^1{I\over(t-\eta)}\ dt+\int_0^1{J\over(t-\theta)}}^{constants}\ dt\right)$$
Each integral is solvable with variable substitution like we did before, we get
$$F_{Area}={-1\over2\delta}\left(G(log(\ldots)-log(\ldots))+H(log(\ldots)\cdots\right)$$
