exponenta event banner

Численное решение уравнений

Символьная математическая Toolbox™ предлагает как числовые, так и символьные решатели уравнений. Сравнение числовых и символьных решателей см. в разделе Выбор числового или символьного решателя. Уравнение или система уравнений может иметь несколько решений. Чтобы найти эти решения численно, используйте функцию vpasolve. Для полиномиальных уравнений: vpasolve возвращает все решения. Для неполиномиальных уравнений vpasolve возвращает первое найденное решение. В этих примерах показано, как использовать vpasolve найти решения как полиномиальных, так и неполиномиальных уравнений, и как получить эти решения с произвольной точностью.

Найти все корни полиномиальной функции

Использовать vpasolve найти все решения функции f (x) = 6x7-2x6 + 3x3-8.

syms f(x)
f(x) = 6*x^7-2*x^6+3*x^3-8;
sol = vpasolve(f)
sol = 

(1.0240240759053702941448316563337-0.88080620051762149639205672298326+0.50434058840127584376331806592405i-0.88080620051762149639205672298326-0.50434058840127584376331806592405i-0.22974795226118163963098570610724+0.96774615576744031073999010695171i-0.22974795226118163963098570610724-0.96774615576744031073999010695171i0.7652087814927846556172932675903+0.83187331431049713218367239317121i0.7652087814927846556172932675903-0.83187331431049713218367239317121i)[vpa('1.0240240759053702941448316563337'); - vpa('0.88080620051762149639205672298326') + vpa('0.50434058840127584376331806592405i'); - vpa('0.88080620051762149639205672298326') - vpa('0.50434058840127584376331806592405i'); - vpa('0.22974795226118163963098570610724') + vpa('0.96774615576744031073999010695171i'); - vpa('0.22974795226118163963098570610724') - vpa('0.96774615576744031073999010695171i'); vpa('0.7652087814927846556172932675903') + vpa('0.83187331431049713218367239317121i'); vpa('0.7652087814927846556172932675903') - vpa('0.83187331431049713218367239317121i')]

vpasolve возвращает семь корней функции, как и ожидалось, поскольку функция является многочленом степени семь.

Поиск нулей неполиномиальной функции с помощью диапазонов поиска и начальных точек

График функции f (x) = e (x/7) cos (2x) показывает периодические нули с увеличением откосов в нулевых точках по мере увеличения x.

syms x
h = fplot(exp(x/7)*cos(2*x),[-2 25]);
grid on

Figure contains an axes. The axes contains an object of type functionline.

Использовать vpasolve чтобы найти ноль функции f. Обратите внимание, что vpasolve возвращает только одно решение неполиномиального уравнения, даже если существует несколько решений. При повторных вызовах, vpasolve возвращает тот же результат.

f = exp(x/7)*cos(2*x);
for k = 1:3
  vpasolve(f,x)
end
ans = -7.0685834705770347865409476123789-vpa('7.0685834705770347865409476123789')
ans = -7.0685834705770347865409476123789-vpa('7.0685834705770347865409476123789')
ans = -7.0685834705770347865409476123789-vpa('7.0685834705770347865409476123789')

Чтобы найти несколько решений, задайте опцию 'Random' кому true. Это делает vpasolve случайным образом выбирайте начальные точки. Сведения об алгоритме, выбирающем случайные начальные точки, см. в разделе Алгоритмы на vpasolve страница.

for k = 1:3
  vpasolve(f,x,'Random',true)
end
ans = -226.98006922186256147892598444194-vpa('226.98006922186256147892598444194')
ans = 98.174770424681038701957605727484vpa('98.174770424681038701957605727484')
ans = 52.621676947629036744249276669932vpa('52.621676947629036744249276669932')

Чтобы найти ноль, близкий к x = 10, установите начальную точку на10.

vpasolve(f,x,10)
ans = 10.210176124166828025003590995658vpa('10.210176124166828025003590995658')

Чтобы найти ноль, близкий к x = 1000, установите начальную точку на1000.

vpasolve(f,x,1000)
ans = 999.8118620049516981407362567287vpa('999.8118620049516981407362567287')

Чтобы найти ноль в 15≤x≤25 диапазона, задайте для диапазона поиска значение [15 25].

vpasolve(f,x,[15 25])
ans = 21.205750411731104359622842837137vpa('21.205750411731104359622842837137')

Поиск нескольких нулей в диапазоне [15 25], звонить нельзя vpasolve многократно, потому что он возвращает один и тот же результат для каждого вызова, как показано ранее. Вместо этого задайте диапазон поиска и задайте 'Random' кому true.

for k = 1:3
  vpasolve(f,x,[15 25],'Random',true)
end
ans = 21.205750411731104359622842837137vpa('21.205750411731104359622842837137')
ans = 21.205750411731104359622842837137vpa('21.205750411731104359622842837137')
ans = 16.493361431346414501928877762217vpa('16.493361431346414501928877762217')

Поскольку 'Random' выбирает начальные точки случайным образом, то же решение может быть найдено при последовательных вызовах.

Найти все нули в указанном диапазоне поиска

Создание функции findzeros систематически находить все нули для f в заданном диапазоне поиска в пределах заданного допуска ошибки. Функция начинается с диапазона входного поиска и вызовов vpasolve чтобы найти ноль. Затем он разбивает диапазон поиска на два вокруг нулевого значения и рекурсивно вызывает себя с новыми диапазонами поиска в качестве входных данных, чтобы найти больше нулей.

Здесь функция описывается в разделе.

Объявите функцию с тремя входами и одним выходом. Первый вход - это функция, второй вход - это диапазон, а дополнительный третий вход позволяет указать ошибку между нулем и более высокой и нижней границами, генерируемыми из него.

function sol = findzeros(f,range,err)

Если дополнительный аргумент для допуска ошибок не указан, findzeros наборы err кому 0.001.

if nargin < 2
  err = 1e-3;
end

Найти ноль в диапазоне поиска с помощью vpasolve.

sol = vpasolve(f,range);

Если vpasolve не находит ноль, выход.

if(isempty(sol))
  return

Если vpasolve находит ноль, разбивает диапазон поиска на два диапазона поиска выше и ниже нуля.

else
  lowLimit = sol-err;
  highLimit = sol+err;

Звонить findzeros с нижним диапазоном поиска. Если findzeros возвращает нули, копирует значения в массив решения и сортирует их.

  temp = findzeros(f,[range(1) lowLimit],1);
  if ~isempty(temp)
    sol = sort([sol temp]);
  end

Звонить findzeros с более высоким диапазоном поиска. Если findzeros возвращает нули, копирует значения в массив решения и сортирует их.

  temp = findzeros(f,[highLimit range(2)],1);
  if ~isempty(temp)
    sol = sort([sol temp]);
  end
  return
end
end

Вся функция findzeros является следующим. Сохранить эту функцию как findzeros.m в текущей папке.

function sol = findzeros(f,range,err)
if nargin < 3
  err = 1e-3;
end
sol = vpasolve(f,range);
if(isempty(sol))
  return
else
  lowLimit = sol-err;
  highLimit = sol+err;
  temp = findzeros(f,[range(1) lowLimit],1);
  if ~isempty(temp)
    sol = sort([sol temp]);
  end
  temp = findzeros(f,[highLimit range(2)],1);
  if ~isempty(temp)
    sol = sort([sol temp]);
  end
  return
end
end

Звонить findzeros с диапазоном поиска [15 25] найти все нули в этом диапазоне для f(x) = exp(x/7)*cos(2*x), в пределах допуска по умолчанию.

syms f(x)
f(x) = exp(x/7)*cos(2*x);
sol = findzeros(f,[15 25])'
sol = 

(16.49336143134641450192887776221718.06415775814131112116019945385719.63495408493620774039152114549721.20575041173110435962284283713722.77654673852600097885416452877624.347343065320897598085486220416)[vpa('16.493361431346414501928877762217'); vpa('18.064157758141311121160199453857'); vpa('19.634954084936207740391521145497'); vpa('21.205750411731104359622842837137'); vpa('22.776546738526000978854164528776'); vpa('24.347343065320897598085486220416')]

Получение решений для произвольной точности

Использовать digits для установки точности решений, возвращаемых vpasolve. По умолчанию vpasolve возвращает решения с точностью до 32 значимых фигур.

f = exp(x/7)*cos(2*x);
vpasolve(f)
ans = -7.0685834705770347865409476123789-vpa('7.0685834705770347865409476123789')

Использовать digits повысить точность до 64 значимых цифр. При изменении digitsубедитесь, что текущее значение сохранено для восстановления.

digitsOld = digits;
digits(64)
vpasolve(f)
ans = -7.068583470577034786540947612378881489443631148593988097193625333-vpa('7.068583470577034786540947612378881489443631148593988097193625333')

Затем измените точность решений на 16 значимых цифр.

digits(16)

Решение многомерных уравнений с помощью диапазонов поиска

Рассмотрим следующую систему уравнений.

z = 10 (cos (x) + cos (y)) z = x + y2-0.1x2yx + y-2.7 = 0

График уравнений для 0≤x≤2.5 и 0≤x≤2.5 показывает, что три поверхности пересекаются в двух точках. Для лучшей визуализации графика используйте view. Для масштабирования значений карты цветов используйте caxis.

syms x y z
eqn1 = z == 10*(cos(x) + cos(y));
eqn2 = z == x+y^2-0.1*x^2*y;
eqn3 = x+y-2.7 == 0;
equations = [eqn1 eqn2 eqn3];
fimplicit3(equations)
axis([0 2.5 0 2.5 -20 10])
title('System of Multivariate Equations')
view(69, 28)
caxis([-15 10])

Figure contains an axes. The axes with title System of Multivariate Equations contains 3 objects of type implicitfunctionsurface.

Использовать vpasolve для поиска точки пересечения поверхностей. Функция vpasolve возвращает структуру. Для доступа к x-, y-, и z-значения решения, индексация в структуру.

sol = vpasolve(equations);
[sol.x sol.y sol.z]
ans = (2.3697477224547980.33025227754520212.293354376823228)[vpa('2.369747722454798'), vpa('0.3302522775452021'), vpa('2.293354376823228')]

Для поиска области пространства решения укажите диапазоны поиска переменных. Если указаны диапазоны 0≤x≤1.5 и 1.5≤y≤2.5, то vpasolve функция выполняет поиск в ограниченной области, показанной на рисунке.

Использовать vpasolve чтобы найти решение для этого диапазона поиска. Чтобы опустить диапазон поиска для z, задайте для третьего диапазона поиска значение [NaN NaN].

vars = [x y z];
range = [0 1.5; 1.5 2.5; NaN NaN];
sol = vpasolve(equations, vars, range);
[sol.x sol.y sol.z]
ans = (0.91062661725633361.7893733827436663.964101572135625)[vpa('0.9106266172563336'), vpa('1.789373382743666'), vpa('3.964101572135625')]

Чтобы найти несколько решений, установите 'Random' опция для true. Это делает vpasolve использовать случайные начальные точки для последовательных прогонов. 'Random' можно использовать в сочетании с диапазонами поиска для создания vpasolve использовать случайные начальные точки в пределах диапазона поиска. Поскольку 'Random' выбирает начальные точки случайным образом, то же решение может быть найдено при последовательных вызовах. Звонить vpasolve многократно, чтобы обеспечить поиск обоих решений.

clear sol
range = [0 3; 0 3; NaN NaN];
for k = 1:5
  temp = vpasolve(equations,vars,range,'Random',true);
  sol(k,1) = temp.x;
  sol(k,2) = temp.y;
  sol(k,3) = temp.z;
end
sol
sol = 

(2.3697477224547980.33025227754520212.2933543768232282.3697477224547980.33025227754520212.2933543768232282.3697477224547980.3302522775452022.2933543768232280.91062661725633361.7893733827436663.9641015721356250.91062661725633361.7893733827436663.964101572135625)[vpa('2.369747722454798'), vpa('0.3302522775452021'), vpa('2.293354376823228'); vpa('2.369747722454798'), vpa('0.3302522775452021'), vpa('2.293354376823228'); vpa('2.369747722454798'), vpa('0.330252277545202'), vpa('2.293354376823228'); vpa('0.9106266172563336'), vpa('1.789373382743666'), vpa('3.964101572135625'); vpa('0.9106266172563336'), vpa('1.789373382743666'), vpa('3.964101572135625')]

Постройте график уравнений. Наложение решений в виде графика рассеяния точек желтым цветом X маркеры с использованием scatter3. Чтобы лучше визуализировать график, сделайте две поверхности прозрачными с помощью alpha. Масштабирование карты цветов до значений графика с помощью caxisи изменить перспективу с помощью view.

vpasolve находит решения на пересечении поверхностей, образованных уравнениями, как показано на рисунке.

clf
ax = axes;
h = fimplicit3(equations);
h(2).FaceAlpha = 0;
h(3).FaceAlpha = 0;
axis([0 2.5 0 2.5 -20 10])
hold on
scatter3(sol(:,1),sol(:,2),sol(:,3),600,'yellow','X','LineWidth',2)
title('Randomly Found Solutions in Specified Search Range')
cz = ax.Children;
caxis([0 20])
view(69,28)
hold off

Figure contains an axes. The axes with title Randomly Found Solutions in Specified Search Range contains 4 objects of type implicitfunctionsurface, scatter.

Наконец, восстановить старое значение digits для дальнейших расчетов.

digits(digitsOld)