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

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

Можно использовать 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')