Документация

Объявите уравнения

Можно использовать Z-преобразование, чтобы решить разностные уравнения, такие как известная "проблема" Роста Кролика. Если пара кроликов назревает за один год, и затем производит другую пару кроликов каждый год, популяция кроликов p (n) в год, n описан этим разностным уравнением.

p (n +2) = p (n +1) + p (n).

Объявите уравнение как выражение, принимающее, что правой стороной является 0. Поскольку n представляет годы, примите, что n является положительным целым числом. Это предположение упрощает результаты.

syms p(n) z
assume(n>=0 & in(n,'integer'))
f = p(n+2) - p(n+1) - p(n)

f =
p(n + 2) - p(n + 1) - p(n)

Решите уравнения

Найдите Z-преобразование уравнения.

fZT = ztrans(f,n,z)

fZT =
z*p(0) - z*ztrans(p(n), n, z) - z*p(1) + z^2*ztrans(p(n), n, z) - z^2*p(0) - ztrans(p(n), n, z)

Функциональный solve решает только для символьных переменных. Поэтому, чтобы использовать solve, сначала замените ztrans(p(n),n,z) с переменными pZT.

syms pZT
fZT = subs(fZT,ztrans(p(n),n,z),pZT)

fZT =
z*p(0) - pZT - z*p(1) - pZT*z - z^2*p(0) + pZT*z^2

Решите для pZT.

pZT = solve(fZT,pZT)

pZT =
-(z*p(1) - z*p(0) + z^2*p(0))/(- z^2 + z + 1)

Вычислите p (n) путем вычисления обратного Z-преобразования pZT. Упростите результат.

pSol = iztrans(pZT,z,n);
pSol = simplify(pSol)

pSol =
2*(-1)^(n/2)*cos(n*(pi/2 + asinh(1/2)*1i))*p(1) + ...
                (2^(2 - n)*5^(1/2)*(5^(1/2) + 1)^(n - 1)*(p(0)/2 - p(1)))/...
                 5 - (2*2^(1 - n)*5^(1/2)*(1 - 5^(1/2))^(n - 1)*(p(0)/2 - p(1)))/5

Замените значениями

Чтобы построить результат, сначала замените значениями начальных условий. Позвольте p(0) и p(1) быть 1 и 2, соответственно.

pSol = subs(pSol,[p(0) p(1)],[1 2])

pSol =
4*(-1)^(n/2)*cos(n*(pi/2 + asinh(1/2)*1i)) - (3*2^(2 - n)*5^(1/2)*(5^(1/2) + 1)^(n - 1))/10 + (3*2^(1 - n)*5^(1/2)*(1 - 5^(1/2))^(n - 1))/5

Постройте результаты

Покажите рост в популяции кроликов в зависимости от времени путем графического вывода pSol.

nValues = 1:10;
pSolValues = subs(pSol,n,nValues);
pSolValues = double(pSolValues);
pSolValues = real(pSolValues);
stem(nValues,pSolValues)
title('Rabbit Population')
xlabel('Years (n)')
ylabel('Population p(n)')
grid on

Анализ результатов

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

Поскольку все функции в pSol могут быть выражены с точки зрения exp, перезапись pSol к exp. Упростите результат при помощи simplify с 80 дополнительные шаги упрощения. Теперь, можно анализировать pSol.

pSol = rewrite(pSol,'exp');
pSol = simplify(pSol,'Steps',80)

pSol =
(2*2^n)/(- 5^(1/2) - 1)^n - (3*5^(1/2)*(1/2 - 5^(1/2)/2)^n)/10 + (3*5^(1/2)*(5^(1/2)/2 + 1/2)^n)/10 - (3*(1/2 - 5^(1/2)/2)^n)/2 + (5^(1/2)/2 + 1/2)^n/2

Визуально осмотрите pSol. Заметьте, что pSol является суммой условий. Каждый термин является отношением, которое может увеличиться или уменьшиться, как n увеличивается. Для каждого термина можно подтвердить эту гипотезу несколькими способами:

Для простоты используйте третий подход. Вычислите условия в n = 100, и затем проверьте подход. Во-первых, найдите отдельные условия при помощи children, замены для n, и преобразуйте, чтобы удвоиться.

pSolTerms = children(pSol);
pSolTermsDbl = subs(pSolTerms,n,100);
pSolTermsDbl = double(pSolTermsDbl)

pSolTermsDbl =
   1.0e+20 *
    0.0000   -0.0000    5.3134   -0.0000    3.9604

Результат показывает, что некоторыми условиями является 0, в то время как другие условия имеют большое значение. Выдвиньте гипотезу, что условия большого значения производят экспоненциальное поведение. Аппроксимируйте pSol с этими условиями.

idx = abs(pSolTermsDbl)>1; % use arbitrary cutoff
pApprox = pSolTerms(idx);
pApprox = sum(pApprox)

pApprox =
(3*5^(1/2)*(5^(1/2)/2 + 1/2)^n)/10 + (5^(1/2)/2 + 1/2)^n/2

Проверьте гипотезу путем графического вывода ошибки приближения между pSol и pApprox. Как ожидалось ошибка переходит к 0, когда n увеличивается. Этот результат демонстрирует, как символьные вычисления помогают вам анализировать свою проблему.

Perror = pSol - pApprox;
nValues = 1:30;
Perror = subs(Perror,n,nValues);
stem(nValues,Perror)
xlabel('Years (n)')
ylabel('Error (pSol - pApprox)')
title('Error in Approximation')