程序的作用是求沙漏重力的变化，因此需要的得到质心的函数，并对质心对时间求二阶导数，即可得到重力变化。
现在质心变化的函数已经编写出来了和曲线也可以画出来（即函数中Z[t]函数，在程序的最后一行）。
但是，第二步对于质心函数对时间求二阶导数部分出现了问题，因为质心的函数是一个分段函数，而且十分复杂，MMA只能给出表达式的解，并不能给出具体的数值，也不能画出曲线图像。
我之后的计划是想将质心曲线上的坐标读出来，然后扔到MATLAB或者别的软件里，通过对离散值求二阶导数 然后得到重力变化。
因此本帖想问的问题是：如何获取曲线上点的坐标？
（PS：曾经尝试用循环语句输出不同的Z[t]值，当时是使t=0,步长为0.1,结果算了两个多小时都没有算出来。）
以下为程序：
rou = 5230; R1 = 0.0275; R2 = 0.01; z1 = 0.082; z2 = 0.035; a =
28.3 \[Degree]; F = 0.69; g = -9.81; R3 = 0.0275; z3 = -0.215;
mi0 = rou*Pi*Power[R1, 3]*Tan[a]/3;
ti0 = 0;
tie = mi0/F;
mie = 0;
zi0 = z1 - R1*Tan[a]/4;
zie = z1 - R1*Tan[a]/3;
mi[t_] := If[t >= ti0, If[t <= tie, mi0 - F*t, mie], mi0];
zi[t_] :=
If[t >= ti0, If[t <= tie, (zi0*mi0 + zx[t]*mx[t])/mi[t], zie], zi0];
zx[t_] := z1 - Power[3*F*t*Power[Tan[a], 2]/(rou*Pi), 1/3]/4;
zxx[t_] = z1 - 3*F*t/(4*rou*Pi*R1^2);
mx[t_] := -F*t;
tii0 = tie;
tiie = tii0 + mii0/F;
mii0 = Pi*Power[R1, 2]*(z1 - z2)*rou;
miie = 0;
zii0 = (z1 + z2)/2 - R1*Tan[a]/3;
ziie = z2 - R1*Tan[a]/3;
mii[t_] :=
If[t >= tii0, If[t <= tiie, mii0 - F*(t - tii0), miie], mii0];
zii[t_] :=
If[t >= tii0, If[t <= tiie, (de2z[t] + z2)/2 - R1*Tan[a]/3, ziie],
zii0];
de2z[t_] := z1 - F*(t - tii0)/(rou*Pi*R1^2);
T = 1/Tan[a] - (R1 - R2)/z2;
m1x = rou*Pi/(3*Tan[a]^2)*(R2/T)^3;
m2x = -rou*Pi/3*(R2/T + z2*R2/(R1 - R2))*(R2/(T*Tan[a]))^2;
m3x = -Pi*rou/3*(R2 + (R1 - R2)/z2*de3z)^3*Tan[a];
m4x = Pi*rou/3*(R2 + (R1 - R2)/z2*de3z)^2*(de3z + z2*R2/(R1 - R2));
z1x = 3*R2/(4*T);
z2x = R2/T - (R2/T + z2*R2/(R1 - R2))/4;
z3x = de3z - Tan[a]*(R2 + (R1 - R2)/z2*de3z)/4;
z4x = 3*de3z/4 - z2*R2/(4*(R1 - R2));
miii0 = m1x + m2x - Pi*rou/3*R1^3*Tan[a] +
Pi*rou/3*R1^2*z2*(1 + R2/(R1 - R2));
miiie = 0;
ziii0 = (z1x*m1x +
z2x*m2x + (z2 - Tan[a]/4*(R2 + (R1 - R2)/z2*z2))*(-Pi*rou*R1^3*
Tan[a]/3) + (3*z2/4 - z2*R2/(4*(R1 - R2)))*(Pi*rou/3*R1^2*
z2*(1 + R2/(R1 - R2))))/miii0;
ziiie = ziii[tiiie];
tiii0 = tiie;
tiiie = tiii0 + miii0/F;
miii[t_] :=
If[t >= tiii0, If[t <= tiiie, miii0 - F*(t - tiii0), miiie], miii0];
ziii[t_] :=
If[t >= tiii0,
If[t <= tiiie, (z1x*m1x + z2x*m2x + z3x*m3x + z4x*m4x)/(miii0 -
F*(t - tiii0)), ziiie], ziii0];
de3z = -0.02 + 5.7405*10^-2 *(4.28891 - 3.6453919 *t)^(1/3);
miv0 = -((R2^3 rou (-R1 + R2 - 2 T z2)*Pi)/(3 T^2 z2));
mive = 0;
tiv0 = tiiie;
tive = tiv0 + miv0/F;
miv[t_] :=
If[t >= tiv0, If[t <= tive, miv0 - F*(t - tiv0), mive], miv0];
ziv0 = 0.0045;
zive = 0;
k = ziv0/(miv0/F);
ziv[t_] :=
If[t >= tiv0, If[t <= tive, ziv0 - k*(t - tiv0), zive], ziv0];
mjet[zl_, zu_] := F/g*(Sqrt[vx^2 + 2*g*zu] - Sqrt[vx^2 + 2*g*zl]);
zle = z3 + R1^2/R3^2*(z1 - z2) + z2/3*(R1^2 + R1*R2 + R2^2)/R3^2;
vx = -F/(rou*Pi*R2^2);
tv0 = 0;
tvi0 = (-Sqrt[vx^2 + 2*g*z3] - vx)/g;
tvie = tve;
tve = tive + (-Sqrt[vx^2 + 2*g*zle] - vx)/g;
zl[t_] :=
2.192851*10^-57 (-9.71052*10^55 - Sqrt[
2.89939*10^109 - 7.260042*10^108 t] + 2.53234*10^55 t);
zux[t_] := vx*t + g*t^2/2;
zlx[t_] := -zux[t]*Pi*R2^2/(Pi*R3^2);
mv[t_] :=
If[t >= tv0,
If[tv0 <= t <= tvi0, mjet[g*t^2/2 + vx*t, 0],
If[tvi0 <= t <= tive, mjet[zl[t], 0],
If[tive <= t <= tve,
mv[tive] - mv[tive]/(tve - tive)*(t - tive)]]], 0];
zjet[zl_, zu_] := (
F ((vx^2 - g zl) Sqrt[vx^2 + 2 g zl] + (-vx^2 + g zu) Sqrt[
vx^2 + 2 g zu]))/(3 g^2*mjet[zl, zu]);
zce = (zle + z3)/2;
zv[t_] :=
If[t >= tv0,
If[tv0 <= t <= tvi0, zjet[g*t^2/2 + vx*t, 0],
If[tvi0 <= t <= tive, zjet[zl[t], 0],
If[tive <= t <= tve,
zv[tive] - (zv[tive] - zce)/(tve - tive)*(t - tive)]]], 0];
mvi0 = 0;
mvie = mi0 + mii0 + miii0 + miv0;
zvi0 = z3;
zvie = (zle + z3)/2;
mvi[t_] :=
If[t >= tvi0, If[t <= tve, rou*Pi*R3^2*(zl[t] - z3), mvie], mvi0];
zvi[t_] := If[t >= tvi0, If[t <= tve, (zl[t] + z3)/2, zvie], zvi0];
M = mi0 + mii0 + miii0 + miv0;
Z[t_] := (mi[t]*zi[t] + mii[t]*zii[t] + miii[t]*ziii[t] +
miv[t]*ziv[t] + mv[t]*zv[t] + mvi[t]*zvi[t])/M;
Plot[Z[t], {t, 0, tve}]

……你看看PiecewiseExpand的帮助。