我有两个二阶 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 的数值解(因为 maple 无法符号求解)。
现在，ODE 在输入中进行了时间变换:
diff(y(t), t, t) = -(diff(y(t), t))-sin(y(t))*cos(y(t))+sin(w*t)
设 B(t) 为该 ODE 的数值解(因为 maple 无法符号解它)。
我的问题是:有没有办法解决不同w值的错误？为此，在对 A(t) 进行数值求解后，我需要将 A(t) 的数值解更改为 A(wt)。我的最终目标是绘制误差与误差的关系图。频率，即 w。当w为1时，应该不会有错误，因为系统没有变化。
我对编码还比较陌生。由于 Maple 无法象征性地解决它，所以我所做的就是以数字方式解决每个问题(使用特定的 w。但是，我想对 [0..1.5] 范围内的 w 进行解决)。然后我将它们绘制在同一个图上。然而，这给了我 A(t) 的数值，而不是 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 之间时误差会增加。
最佳答案
有多种方法可以做到这一点。有些比其他的更有效率。有些更方便。我会展示一些。
第一种基本方法涉及对 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
关于numeric - 在 ODE 求解后，如何将数值解更改为 ODE？，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/56451293/