vpasolve (Matlab) не может найти root уравнения при кусочном использовании (root явно существует)

Question

Я пытаюсь найти root этого уравнения, которое является только функцией H:

syms H
eq = piecewise(H < 0.624, - 98.3*H^6 - 0.961*H^3 + 6.94*H^2 + 2.41e-8, H < 1.56, 6.94*H^2 - 5377.0*(0.641*H - 0.4)^3 - 5077.0*(0.641*H - 0.4)^2 - 0.961*H^3 - 98.3*H^6 + 2.41e-8, 1.56 <= H, - 98.3*H^6 - 0.961*H^3 + 6.94*H^2 - 2988.0);
H = vpasolve(eq==0,H);

Однако vpasolve возвращает пустой символ 0x1. Если я построю график с помощью:

fplot(eq2), xlim([0.5 0.51]), ylim([-0.1 +0.1])

, вы ясно увидите, что он пересекает ось y = 0 около 0,506. введите описание изображения здесь Фактически:

double(subs(eq,0.5062))

Дает -8.2450e-05.

виновата, кажется, кусочная функция. Ясно, что первое условие удовлетворяется H, и решение является решением - 98,3 H ^ 6 - 0,961 H ^ 3 + 6,94 * H ^ 2 + 2,41e-8. Если удалить кусочно, решатель работает:

syms H
eq - 98.3*H^6 - 0.961*H^3 + 6.94*H^2 + 2.41e-8;
H = vpasolve(eq==0,H);

, что дает

H =

                                                 -0.52458969841163012770887662427408
                                                  0.50619380851052153463867425194323
- 0.00000000024038304755545478229904837943596 + 0.000058924789608363053116092886242865i
- 0.00000000024038304755545478229904837943596 - 0.000058924789608363053116092886242865i
          0.0091979451909373440905559684644729 - 0.51555587519744206224638650619188i
          0.0091979451909373440905559684644729 + 0.51555587519744206224638650619188i

Положительное действительное решение является искомым. А вообще мне нужна кусочная, почему не работает решатель? Обходные пути?

Изменить:

Некоторые из вас предлагают использовать «решение», которое действительно работает. Я действительно пробовал решить перед публикацией, но я исключил его, потому что он в 5 раз медленнее, чем vpasolve (тестировался только на первом уравнении, поэтому без кусочков). Мне нужно решать это тысячи раз для разных параметров. было бы замечательно, если бы этот баг с vpasolve можно было решить, нигде не сказано, что vpasolve не принимает кусочные функции

Answer 1

Я не уверен, что vpasolve совместим с кусочной функцией, но я уверен, что solve совместим, поэтому вы можете использовать:

res = vpa(solve(eq==0,H))

И мы получаем два корня:

res =  0.5026 
      -0.5245

Вы уверены, что установили правильное условие?

H < 0.624 -> eq1
H < 1.56  -> eq2
H >= 1.56 -> eq3

eq1 и eq2 перекрываются, я бы использовал 0.624 <= H < 1.56 -> eq2, чтобы прояснить ситуацию. И если вам нужно только положительное и реальное решение, вы также должны добавить это условие прямо в свое уравнение:

0     <= H < 0.624  -> eq1
0.624 <= H < 1.560  -> eq2
H     >=     1.560  -> eq3

И использовать assume:

assume(H,'real')

Итак, теперь все идеально ясно.