2017-01-09 10 views
1

Я пытаюсь узнать, как использовать MATLAB для решения системы дифференциальных уравнений (уравнения Лоренца) и сюжет каждого решения в зависимости от тРешение системы ОДУ с помощью ode45

X’ = −σx + σy 
Y’ = ρx − y − xz 
Z’ = −βz + xy 

где σ = 10, β = 8/3 , и ρ = 28, а также x(0) = −8, y(0) = 8 и z(0) = 27.

Вот код, который я использую:

function xprime = example(t,x) 

sig = 10; 
beta = 8/3; 
rho = 28; 
xprime = [-sig*x(1) + sig*x(2); 
      rho*x(1) - x(2) - x(1)*x(3); 
      -beta*x(3) + x(1)*x(2)]; 

x0 = [-8 8 27];  
tspan = [0 20];  
[t,x] = ode45(@example, tspan, x0); 

figure  
plot(t,x(:,1)), hold on 
plot(t,x(:,2)), hold on 
plot(t,x(:,3)), hold off 

Однако это порождает ошибку, как это исправить? Я не уверен, какие входные аргументы отсутствуют или я ошибаюсь. Я ценю любую помощь, спасибо.

Недостаточно входных аргументов.

Ошибка в примере (строка 9) xprime = [- sig x (1) + sig x (2); rho * x (1) - x (2) - x (1) x (3); -beta x (3) +
x (1) * x (2)];

+0

Вы рекурсии - это то, что намеренно, или ...? –

+0

Возможно, я только это, потому что я не знаю, как лучше это написать. – Tina

ответ

3

Это на самом деле довольно хорошо с первой попытки!

Проблема заключается в том, что при нажатии кнопки Run (или нажмите F5) вы вызываете функцию example без аргументов; о чем жалуется MATLAB.

Вторая проблема заключается в том, что, даже если вы должны были быть в состоянии выполнить функцию, как это, ode45 бы вызвать функцию example, которая будет вызывать ode45, который назвал бы example, который назвал бы ode45 и так далее, до тех пор, предел рекурсии.

Решение как это разделить его на две функции (они могут быть записаны в тот же М-файл):

% top-level function; no arguments 
function caller() 

    x0 = [-8 8 27]; 
    tspan = [0 20]; 
    [t,x] = ode45(@example, tspan, x0); 

    figure 
    plot(t,x(:,1)), hold on 
    plot(t,x(:,2)), hold on 
    plot(t,x(:,3)), hold off 

end 

% The derivative 
function xprime = example(t,x) 

    sig = 10; 
    beta = 8/3; 
    rho = 28; 
    xprime = [-sig*x(1) + sig*x(2); 
       rho*x(1) - x(2) - x(1)*x(3); 
       -beta*x(3) + x(1)*x(2)];    
end