Существует множество способов сделать это. Некоторые из них более эффективны, чем другие. Некоторые из них более удобны. Я покажу несколько.
Этот первый основной способ включает два вызова 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 как параметра, который мы настраиваем на лету. (Мы не хотим нести расходы или неудобства при вызове dsolve для каждого значения w.)
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 текущему сохраненному значению параметра и, если он отличается, устанавливает хранимое значение. Затем он вызывает вычислительную процедуру Ygen с заданным значением t.
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]);
Мы также можем сделать (соединенные линии) в точечных графиках, что более эффективно, поскольку все значения t вычисляются последовательно для любого заданного значения w. (То есть, оно не колеблется между 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) Используйте специальную процедуру Yw1 для w = 1.0, чтобы Yw не вызывался для установки параметра w дважды для каждой оценки подынтегрального выражения (в норме).
2) Используйте опцию запомнить эту процедуру Yw1.
3) Используйте опцию compile=true в двух dsolve вызовах.
Я также исправляю формулу в норме для вызова 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
