Я предполагаю, что под IVF(V1(t)) вы действительно имели в виду F(V(t)).
Мне показалось, что вычисление ваших dV(t) (и, возможно, V(t) иногда) для очень маленьких tможет иметь некоторые числовые трудности.
И я волновался, что использование более узкого диапазона интеграции, например, скажем, t=1.0e-2 .. 1.0, может привести к некоторому значительному вкладу в интеграл из этой части.
Итак, ниже я построил V и dV так, чтобы они повысили Digits (только) для большей точности обработки (только) собственного тела вычислений.
И я переключил его ссимволическая интеграция в числовую / поплавковую интеграцию.И для этого я добавил опцию epsilon на тот случай, если вам нужно более точно контролировать допуск на точность.
Некоторые из этих деталей помогают решить проблемы со скоростью, о которых вы упоминали выше.И некоторые из этих деталей немного замедляют вычисления, хотя, возможно, и дают большее ощущение их надежности.
Приведенные ниже вызовы команды plotting приведены только для иллюстрации.Прокомментируйте их, если хотите.
Я использовал Maple 2018.0 для них.Если у вас возникли трудности, сообщите о вашей конкретной версии.
restart;
T:=sqrt(x)/(sqrt(x)+sqrt(1-x))^2:
Vexpr:=solve(t=T,x,explicit)[1]:
Vexpr:=simplify(combine(Vexpr)) assuming t>0, t<1:
Процедура вычисления V построена следующим образом.Он построен таким тонким способом по двум причинам:
1) Он внутренне повышает рабочую точность Digits, чтобы избежать ошибки округления.
2) Создает (и выбрасывает) пустойсписок [].В результате он будет не работать в режиме быстрой аппаратной точности evalhf.Это ключ позже, во время числовой интеграции.
V:=subs(__dummy=Vexpr,
proc(t)
if not type(t,numeric) then
return 'procname'(t);
end if;
[];
Digits:=100;
evalf(__dummy);
end proc):
V(0.25);
0.2037579200498004992002294012453548811286405373\
653346413644953624868320875151070347969077227713469370
V(1e-9);
1.0000000040000000119998606410199120814521886485524\
-18
22617659574502917213370846924673012122642160567565 10
Давайте построим T и V,
plot([T, V(x)], x=0..1);
Теперь давайте создадим процедуру для производной V, используя ту же технику, чтобы она тоже подняла Digits (только для себя) и не работала с evalhf.
dV:=subs(__dummy=diff(Vexpr,t),
proc(t)
if not type(t,numeric) then
return 'procname'(t);
end if;
[];
Digits:=100;
evalf(__dummy);
end proc):
Обратите внимание, что вызов dV(t) для неизвестного символа t возвращает его без оценки.Это удобно позже.
dV(t);
dV(t)
dV(0.25);
2.44017580084567947538393626436824494366329948208270464559139762\
2347525580165201957710520046760103982
dV(1e-15);
-5.1961525404198771358909147606209930290335590838862019038834313\
73611362758951860627325613490378754702
evalhf(dV(1e-15)); # I want this to error-out.
Error, unable to evaluate expression to hardware floats: []
Это ваш F.
F:=t->-10*ln(1-t*(1-exp(-1))):
Теперь создайте процедуру U.Это также возвращает неоцененное значение, когда любой из его аргументов не является действительным числом.Обратите внимание, что lambda является его вторым параметром.
U:=proc(t,lambda)
if not ( type(t,numeric) and
type(lambda,numeric) ) then
return 'procname'(args);
end if;
piecewise(t <= .17215,
min(F(V(t)), max(0, 12-(1/4)/(lambda^2*dV(t)^2))),
t >= .17215,
min(F(V(t)), max(0, 12-(1/4)*.7865291304*lambda^(-2))));
end proc:
Давайте попробуем определенное значение t=0.2 и lambda=2.5.
evalf(U(0.2, 2.5));
.6774805135
Давайте построим его.Это займет немного времени.Обратите внимание, что для небольших t.
plot3d([Re,Im](U(t,lambda)),
t=0.0001 .. 1.0, lambda=-2.0 .. 2.0,
color=[blue,red]);
кажется, что он не взорвется. Давайте также построим график dV.Это занимает некоторое время, так как dV повышает Digits.Обратите внимание, что он выглядит точным (не взрывается) для небольших t.
plot(dV, 0.0 .. 1.0, color=[blue,red], view=-1..4);
Теперь давайте создадим процедуру H, которая вычисляетчисловой интеграл.
Обратите внимание, что lambda - это его первый параметр, который он передает в качестве второго аргумента U.Это относится к вашему последующему комментарию.
Второй параметр H передается evalf(Int(...)) для определения его погрешности точности.
Обратите внимание, что диапазон интегрирования составляет от 0,001 до 1,0 (дляt).Вы можете попробовать более низкую конечную точку ближе к 0,0, но это займет еще больше времени для вычислений (и может не увенчаться успехом).
Существует сложное взаимодействие между нижней конечной точкой числовой интеграции и допуском точностиeps и числовое поведение для U и V и dV для очень маленьких t.
(Рабочая точность Digits установлена в 15 внутри H.просто чтобы можно было использовать метод _d01ajc. Тем не менее, V и dV будут повышать точность работы для своих собственных вычислений при каждом числовом вызове подынтегрального выражения. Также я использую unapplyи принудительный выбор метода числовой интеграции в качестве хитрого способа не допустить, чтобы evalf(Int(...)) тратил впустую много времени, пытаясь найти сингулярность, чтобы найти особенности. Вам не нужно полностью понимать все эти детали. Я просто объясняю, что ябыли причины для этой сложной установки.)
H:=proc(lambda,eps)
local res,t;
Digits:=15;
res:=evalf(Int(unapply(U(t,lambda)*dV(t),t),
1e-3..1,
epsilon=eps, method=_d01ajc) - R);
if type(res,numeric) then
res;
else
Float(undefined);
end if;
end proc:
Давайте попробуем позвонить H для конкретного R.Обратите внимание, что это не быстро.
R:=2.93;
R := 2.93
H(0.1,1e-3);
-2.93
H(0.2,1e-5);
0.97247264305333
Давайте построим график H для фиксированного допуска интеграции.Это занимает довольно много времени.
[отредактировано] Я изначально подготовил t->H(t,1e-3) здесь.Но ваш последующий комментарий показал недопонимание цели первого параметра этого оператора.t в t->H(t,1e-3) является просто фиктивным именем и НЕ относится к параметру t V и dV, второму параметру U и интегрируемости.Важно то, что этот фиктивный объект передается в качестве первого аргумента H, который, в свою очередь, передает его в качестве второго аргумента U, который используется для lambda.Поскольку раньше была путаница, сейчас я буду использовать фиктивное имя LAMBDA.
plot(LAMBDA->H(LAMBDA,1e-3), 0.13 .. 0.17 , adaptive=false, numpoints=15);
Теперь давайте попробуем найти значениеэтот корень (перехват).
[отредактировано] Здесь я также ранее использовал фиктивное имя t для первого аргумента H.Это просто фиктивное имя, которое не относится к параметру t V и dV, второму параметру U и переменной интеграции.Поэтому я изменю его ниже на LAMBDA, чтобы уменьшить путаницу.
Это не быстро.Я минимизирую абсолютное значение, потому что проблемы с числовой точностью делают это проблематичным для числового корневого поиска fsolve.
Osol:=Optimization:-Minimize('abs(H(LAMBDA,1e-3))', LAMBDA=0.15..0.17);
[ -8 ]
Osol := [1.43561000000000 10 , [LAMBDA = 0.157846355888300]]
Мы можем выбрать числовое значение лямбды (которое «решило» интеграцию).
rhs(Osol[2,1]);
0.157846355888300
Если мы вызовем H, используя это значение, то получимостаток (аналогично первой записи в Osol, но без учета абсолютного значения).
H(rhs(Osol[2,1]), 1e-3);
-8
-1.435610 10
Теперь у вас есть инструменты для более точного получения корня.Для этого вы можете поэкспериментировать с:
- жестче (меньше)
eps допуск численного интегрирования. - Нижняя конечная точка интегрирования ближе к 0,0
- Еще выше
Digits принудительно в процедурах V и dV.
Это потребует терпения, поскольку вычисления занимают некоторое время, а числовая интеграция может дать сбой (принимаянавсегда) когда эти три аспекта не работают должным образом вместе.
Я отредактировал тело процедуры H и изменил диапазон интегрирования на 1e-4 .. 1.0, тогда я получаю следующий результат (указав 1e-5 для eps в вызове H).
Osol:=Optimization:-Minimize('abs(H(LAMBDA,1e-5))', LAMBDA=0.15..0.17);
[ -7 ]
Osol := [1.31500690000000 10 , [LAMBDA = 0.157842787382700]]
Таким образом, вполне вероятно, что эта приблизительная корневая лямбда точна с точностью до 4 или 5 десятичных знаков.
Удачи.