如何绘制特殊点的实验数据?

How to plot experimental data at special points?

我有一组实验数据(P),想获取plot"experimental vs predicted"。为此,我使用另一组依赖于 P (Q) 的数据,绘制散点图,使用适当的拟合,然后获得回归线,并在适当的微分方程中使用其系数。 P 的情节看起来不错,但我需要在那里添加实验数据。为简单起见,我使用了区间 t=0..150

如何绘制实验数据以便 P(0) = Pvals[1], P(10)=Pvals[2] 等?此外,我如何分发数据(比如,我有 t=0..800 并且想要绘制 Pvals 以便 P(0) = Pvals[1] and P(800) = Pvals[16])?

Pvals := [3.929, 5.308, 7.24, 9.638, 12.866, 17.069, 23.192, 31.433, 38.558, 50.156, 62.948, 
75.996, 91.972, 105.711, 122.775, 131.669]:

for i to 15 do Qval[i] := .1*(Pvals[i+1]/Pvals[i]-1); end do:
Qvals := [seq(Qval[i], i = 1 .. 15), 0.144513895e-1]:               
    with(Statistics);
ScatterPlot(Pvals, Qvals, fit = [a*v^2+b*v+c, v], thickness = 3, 
legend = [points = "Point data", fit = typeset("fit to a", 2^nd, "degree polynomial")]);

with(CurveFitting);
LeastSquares(Pvals, Qvals, v, curve = a*v^2+b*v+c);

de := diff(P(t), t) = (0.370152282598477e-1-0.272504103112702e-3*P(t))*P(t);

sol := dsolve({de, P(0) = 3.929}, P(t));

P := plot(rhs(sol), t = 0 .. 160);

我不确定我是否完全遵循您的方法。但这是否像您想要完成的那样?

restart;
with(Statistics):

Pvals := [3.929, 5.308, 7.24, 9.638, 12.866, 17.069, 23.192, 31.433,
          38.558, 50.156, 62.948, 75.996, 91.972, 105.711, 122.775, 131.669]:

for i to 15 do Qval[i] := .1*(Pvals[i+1]/Pvals[i]-1); end do:
Qvals := [seq(Qval[i], i = 1 .. 15), 0.144513895e-1]:

form := a*v^2+b*v+c:

CF := CurveFitting:-LeastSquares(Pvals, Qvals, v, curve = form);

       CF := 0.0370152282598477 - 0.000272504103112702 v

                               -7  2
              + 5.60958249026713 10   v

现在我在 DE 中使用 CF(因为我不明白你为什么去掉 v^2 项),

#de := diff(P(t), t) = (0.370152282598477e-1-0.272504103112702e-3*P(t))*P(t);
de := diff(P(t), t) = eval(CF, v=P(t))*P(t);

           d         /                                              
    de := --- P(t) = [=11=].0370152282598477 - 0.000272504103112702 P(t)
           dt                                                       

                            -7     2\     
       + 5.60958249026713 10   P(t) / P(t)

我将使用 dsolve 命令的 numeric 选项,并获得一个计算 P(t) 数值 t 值的过程。

sol := dsolve({de, P(0) = 3.929}, P(t), numeric, output=listprocedure ):

Pfunc := eval(P(t), sol);

              Pfunc := proc(t)  ...  end;

Pfunc(0.0), Pvals[1];

                3.92900000000000, 3.929

现在进行一些重新调整(这也是我对您目标的猜测),

endpt := fsolve(Pfunc(t)-Pvals[16]);

                  endpt := 135.2246055

Pfunc(endpt), Pvals[16];

               131.669000003321, 131.669

plot(Pfunc(t), t=0 .. endpt, size=[500,200]);

a,b,N := 0.0, 800.0, nops(Pvals);

                a, b, N := 0., 800.0, 16

Pfuncscaled := proc(t) 
                 if not t::numeric then
                   return 'procname'(args);
                 end if;
                 Pfunc(t*endpt/b);
               end proc:

Pfuncscaled(0), Pvals[1];

                3.92900000000000, 3.929

Pfuncscaled(800), Pvals[N];

               131.669000003321, 131.669

PLscaled := plot( Pfuncscaled(t), t=a .. b,
                  color=red, size=[500,200] );

现在还要针对 0 .. 800 显示 Pdata

V := Vector(N, (i)->a+(i-1)*(b-a)/(N-1)):

V[1], V[-1];

                    0., 800.0000000

Pdatascaled := plot( < V | Vector(Pvals) >,
                     color=blue, size=[500,200],
                     style=pointline, symbol=solidcircle );

并且,显示重新缩放的数据以及来自 dsolve

的重新缩放的过程
plots:-display( PLscaled, Pdatascaled, size=[500,500] );