求解后如何更改 ODE 的数值解？
How to make a change in a numeric solution to an ODE after it has been solved?
我有两个二阶 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)。我的最终目标是绘制错误 v.s。频率，即 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 之间时误差会增加。
有多种方法可以做到这一点。有些比其他的更有效率。有些更方便。我会展示一些。
第一种基本方法涉及两次调用 dsolve，类似于您的原始代码，但带有额外的 output=listprocedure 选项。这使得 dsolve return 成为 scalar-valued 过程的列表，而不是单个过程（return 是值列表）。这允许我们选择 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的过程。在不使用数字参数调用时，它 return 未计算。当使用数字参数调用时，它会检查 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]);
我们也可以在 point-plots 中执行（连接线），这样效率更高，因为对于任何给定的 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) 在该过程中使用选项 remember 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
我有两个二阶 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)。我的最终目标是绘制错误 v.s。频率，即 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 之间时误差会增加。
有多种方法可以做到这一点。有些比其他的更有效率。有些更方便。我会展示一些。
第一种基本方法涉及两次调用 dsolve，类似于您的原始代码，但带有额外的 output=listprocedure 选项。这使得 dsolve return 成为 scalar-valued 过程的列表，而不是单个过程（return 是值列表）。这允许我们选择 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的过程。在不使用数字参数调用时，它 return 未计算。当使用数字参数调用时，它会检查 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]);
我们也可以在 point-plots 中执行（连接线），这样效率更高，因为对于任何给定的 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) 在该过程中使用选项 remember 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