求解后如何更改 ODE 的数值解？

Question

我有两个二阶 ODE，我想使用（使用 Maple）计算它们之间的误差

|| B(t) - A(w*t) || = sqrt(int(B(t) - A(w*t) )^2), t = 0..T)

其中 A(t) 是对输入没有任何时间变换的系统的解，B(t) 是对输入进行时间变换的系统的解（时间变换为 wt（w 通常很小) 并且 T 是一个数字，而不是变量。我会更改 T 以研究不同时间跨度的误差）。

例如（帮助解释我的问题）：

原来的ODE是：

diff(g(t), t, t) = -(diff(g(t), t))-sin(g(t))*cos(g(t))+sin(t)

令 A(t) 为原始 ODE 的 NUMERIC 解（因为 maple 无法符号解）。

现在，输入中有时间变换的 ODE：

diff(y(t), t, t) = -(diff(y(t), t))-sin(y(t))*cos(y(t))+sin(w*t)

令 B(t) 为该 ODE 的 NUMERIC 解（因为 maple 无法以符号方式求解）。

我的问题是：有什么方法可以解决不同 w 值的错误？为此，在数值求解 A(t) 后，我需要将 A(t) 的数值解更改为 A(wt)。我的最终目标是绘制误差与频率，即 w。当 w 为 1 时，应该没有错误，因为系统没有变化。

我对编码还是比较陌生。我所做的，因为 Maple 无法象征性地解决它，而是以数字方式解决每个问题（使用特定的 w。但是，我希望在 [0..1.5] 范围内对 w 执行此操作）。然后我将它们绘制在同一个图上。但是，这给了我 A(t) NOT A(wt) 的数值。而且，我不知道如何减去它们。

sol1 := dsolve([diff(g(t), t, t) = -(diff(g(t), t))- sin(g(t))*cos(g(t))+sin(t), g(0) = 0, (D(g))(0) = 0], numeric);

sol2 := dsolve([diff(y(t), t, t) = -(diff(y(t), t))-sin(y(t))*cos(y(t))+sin(.5*t), y(0) = 0, (D(y))(0) = 0], numeric);

S1 := plots[odeplot](sol1, 0 .. 10, color = red);
S2 := plots[odeplot](sol2, 0 .. 10, color = blue);
display(S1, S2);

但是，这没有帮助，因为它只绘制 A(t)，而不是 A(wt)。同样，它只绘制它，它没有向我显示它们之间的错误。

我预计误差会随着频率 w 接近 0 接近 0。我确实希望当 w 在 0 和 1 之间时误差会增加。

Answer 1

有多种方法可以做到这一点。有些比其他的更有效。有些更方便。我会展示一些。

第一种基本方法涉及对dsolve 的两次调用，类似于您的原始代码，但带有额外的output=listprocedure 选项。这使得dsolve 返回一个标量值过程列表，而不是单个过程（它将返回一个值列表）。这使我们可以选择y(t) 和g(t) 的各个过程并分别使用它们。

restart;
sol1 := dsolve([diff(g(t), t, t)
                = -(diff(g(t), t))- sin(g(t))*cos(g(t))+sin(t),
                g(0) = 0, (D(g))(0) = 0], numeric,
                output=listprocedure):
sol2 := dsolve([diff(y(t), t, t)
                = -(diff(y(t), t))-sin(y(t))*cos(y(t))+sin(.5*t),
                y(0) = 0, (D(y))(0) = 0], numeric,
                output=listprocedure):

如果你愿意，你仍然可以在这里使用plots:-odeplot。

S1 := plots:-odeplot(sol1, 0 .. 20, color = red, numpoints=200):
S2 := plots:-odeplot(sol2, 0 .. 20, color = blue, numpoints=200):
plots:-display(S1, S2, size=[400,150]);

但您也可以提取单个过程，绘制它们，或绘制它们的差异。

G := eval(g(t),sol1):
Y := eval(y(t),sol2):

plot([G,Y], 0..20, color=[red,blue], size=[400,150]);

它们之间的区别现在更容易绘制出来了。

plot(G-Y, 0..20, color=black, size=[400,150]);

我们可以计算一个范数（你的积分），

sqrt(evalf(Int( t -> ( G(t) - Y(t) )^2, 0..20, epsilon=1e-5)));
         8.74648880831735

但现在让我们更方便地将w 视为我们动态调整的参数。 （我们不想为每个w 值调用dsolve 带来成本或不便。）

solgen := dsolve([diff(y(t), t, t)
                 = -(diff(y(t), t))-sin(y(t))*cos(y(t))+sin(w*t),
                 y(0) = 0, (D(y))(0) = 0], numeric,
                 parameters = [w],
                 output=listprocedure):

 Ygen := eval(y(t),solgen):

 # set the parameter value
 Ygen(parameters=[w=0.5]);
                       [w = 0.5]

 # we could query the parameter value
 Ygen(parameters);
                       [w = 0.5]

 # Now call it at a particular value of t
 Ygen(3.2);
                   0.864744459594622

现在我们构造一个w 和t 的过程。当没有数字参数调用时，它返回未评估。当使用数字参数调用时，它会检查 w 值是否与当前存储的参数值匹配，如果不同，则设置存储的值。然后它以给定的t 值调用计算过程Ygen。

Yw := proc(t,w)
  if not [t,w]::list(numeric) then
    return 'procname'(args);
  end if;
  if w - eval(':-w',Ygen(parameters)) <> 0.0 then
    Ygen(':-parameters'=[':-w'=w]);
  end if;
  Ygen(t);
end proc:

这可以产生与以前相同的图。

plots:-display(
  plot(Yw(t,0.5), t=0..20, color=red),
  plot(Yw(t,1.0), t=0..20, color=blue),
  size=[400,150]
);

# somewhat inefficient since for each t value it
# switches the value of the w parameter.
plot(Yw(t,1.0)-Yw(t,0.5), t=0..20, size=[400,150]);

我们还可以在点图中进行（连接线），因为对于任何给定的w 值，所有t 值都是按顺序计算的，因此效率更高。 （即，它不会在 w 值之间反弹。）

a,b := 0,20:
xpts := Vector(100,(i)->a+(b-a)*(i-1)/(100-1),datatype=float[8]):
plots:-display(
  plot(<xpts | map(t->Yw(t,0.5), xpts)>, color=red),
  plot(<xpts | map(t->Yw(t,1.0), xpts)>, color=blue),
  size=[400,150]
);

# more efficient since all the Yw values for w=1.0 are
# computed together, and all the Yw values for w=0.5 are
# computed together.
plot(<xpts | map(t->Yw(t,1.0), xpts) - map(t->Yw(t,0.5), xpts)>,
     size=[400,150]);

evalf([seq( ['w'=w, sqrt(Int( unapply( (Yw(t,1.0)
                                         - Yw(t,w))^2, t),
                              0..20, epsilon=1e-1))],
             w=0..1.0, 0.1 )]);

    [[w = 0., 2.97123678486962], [w = 0.1, 20.3172660599286], 

     [w = 0.2, 21.8005838723429], [w = 0.3, 13.0097728518328], 

     [w = 0.4, 9.28961600039575], [w = 0.5, 8.74643983270251], 

     [w = 0.6, 6.27082651159143], [w = 0.7, 5.38965679479886], 

     [w = 0.8, 5.21680809275065], [w = 0.9, 3.19786559349464], 

     [w = 1.0, 0.]]


plot(w->sqrt(Int( (Yw(t,1.0) - Yw(t,w))^2, t=0..20,
                  epsilon=1e-3, method=_d01ajc )),
     0..1, size=[400,150]);

plot3d( Yw(t,w), w=0..1.0, t=0..50, grid=[179,179] );

# For 3D plotting it's also faster to compute all t-values
# for each given w-value, rather than the other way round.
CodeTools:-Usage(
  plot3d( Yw(t,w), w=0..1.0, t=0..50, grid=[49,49] )
):

   memory used=37.51MiB, alloc change=0 bytes, cpu time=504.00ms,
   real time=506.00ms, gc time=127.31ms

CodeTools:-Usage(
  plot3d( Yw(t,w), t=0..50, w=0..1.0, grid=[49,49] ) ):

   memory used=124.13MiB, alloc change=0 bytes, cpu time=2.66s,
   real time=2.66s, gc time=285.47ms

[编辑] 该规范的上述情节并不快。下面我做了三个调整以获得更好的性能： 1) 对 w=1.0 使用专用过程 Yw1，这样就不会调用 Yw 为每次评估被积函数（在范数中）设置参数 w 两次。 2）使用选项记住该程序Yw1。 3) 在两个dsolve 调用中使用选项compile=true。

我还将规范中的公式更正为调用Yw1(w*t)，以匹配有问题的原始公式。

restart;
solw1 := dsolve([diff(y(t), t, t)
             = -(diff(y(t), t))-sin(y(t))*cos(y(t))+sin(t),
             y(0) = 0, (D(y))(0) = 0], numeric,
             compile=true,
             output=listprocedure):
Yw1raw := eval(y(t),solw1):
Yw1 := proc(t)
  option remember;
  if not t::numeric then return 'procname'(args); end if;
  []; # evalhf error that plot will catch
  Yw1raw(t);
end proc:
solgen := dsolve([diff(y(t), t, t)
             = -(diff(y(t), t))-sin(y(t))*cos(y(t))+sin(w*t),
             y(0) = 0, (D(y))(0) = 0], numeric,
             parameters = [w],
             compile=true,
             output=listprocedure):
Ygen := eval(y(t),solgen):
Yw := proc(t,w)
  if not [t,w]::list(numeric) then
    return 'procname'(args);
  end if;
  if w - eval(':-w',Ygen(parameters)) <> 0.0 then
    Ygen(':-parameters'=[':-w'=w]);
  end if;
  Ygen(t);
end proc:

CodeTools:-Usage(
  plot(w->sqrt(Int( (Yw1(w*t) - Yw(t,w))^2, t=0..20,
                    epsilon=1e-3, method=_d01ajc )),
       0..1, size=[400,150])
);

    memory used=0.76GiB, alloc change=205.00MiB,
    cpu time=8.22s, real time=7.78s, gc time=761.80ms