Есть ли способ построить функцию f с уравнениями (из системы ОДУ) внутри нее? - PullRequest
Новая
       22

Есть ли способ построить функцию f с уравнениями (из системы ОДУ) внутри нее?

1 голос
/ 06 мая 2019

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

Учитывая эту систему ODE: 

dy(1) = y(2) − y(3);
dy(2) = −3*y(2) − 5*sin(ω*y(1));
dy(3) = y(1)*y(2);

с omega=3. Начальные значения: y(1)=0; y(2)=4; y(3)=1;

Участок f между t=0 и t=20: 

f(t) = |2*cos(y1(t)) + y2(t)|

Сначала я попытался смоделировать систему уравнений с ODE45 с помощью этого кода:

В редакторе: 

function dy=modelo10(t,y)
    global omega

    dy = zeros(3,1);
    dy(1)= y(2) - y(3);
    dy(2)= -3*y(2) - 5*sin(omega*y(1));
    dy(3)= y(1)*y(2);

end

На вкладке команд: 

>>global omega;
>>omega = 3;
>>[T,Y] = ode45(@modelo10, [0,20], [0,4,1]); 
%% I assign the function to variable f
>>f = @(t) 2.*cos(y(:,1)) + y(:,2);
%% And finally plot it with the values t
>>fplot(f, [0,20]);

Я получаю пустой график со следующей ошибкой: 

Warning: Function behaves unexpectedly on array inputs. To improve 
performance, properly
vectorize your function to return an output with the same size and 
shape as the input
arguments. 
  In matlab.graphics.function.FunctionLine>getFunction
  In matlab.graphics.function.FunctionLine/updateFunction
  In matlab.graphics.function.FunctionLine/set.Function_I
  In matlab.graphics.function.FunctionLine/set.Function
  In matlab.graphics.function.FunctionLine
  In fplot>singleFplot (line 241)
  In fplot>@(f)singleFplot(cax,{f},limits,extraOpts,args) (line 196)
  In fplot>vectorizeFplot (line 196)
  In fplot (line 166) 
Warning: Error updating FunctionLine.

The following error was reported evaluating the function in 
FunctionLine update:
Non-scalar in Uniform output, at index 1, output 1.
Set 'UniformOutput' to false.

Warning: Error updating FunctionLine.

The following error was reported evaluating the function in 
FunctionLine update:
Non-scalar in Uniform output, at index 1, output 1.
Set 'UniformOutput' to false

Так что мой вопрос в том, каковы шаги (или команды) для построения f?

Я использую MATLAB R2019a.

Ответы [ 2 ]

4 голосов
/ 06 мая 2019

Ваша первая проблема в определении f: 

[T,Y] = ode45(@modelo10, [0,20], [0,4,1]); 
f = @(t) 2.*cos(y(:,1)) + y(:,2);

Обратите внимание, что в первой строке указана переменная Y, но во второй вы используете y. Это разные переменные. Это, вероятно, работает для вас, потому что у вас есть переменная y, определенная в вашем рабочем пространстве, но она делает не то.

Далее ваш f является функцией f(t), но он определен без использования ввода t. Он всегда выводит одно и то же. fplot ожидает функцию, которая выводит массив того же размера, что и его вход, это сообщение об ошибке, которое вы получаете.

Однако вам не нужно здесь определять функцию, вы можете напрямую вычислить все значения для f и построить их, используя plot: 

[t,y] = ode45(@modelo10, [0,20], [0,4,1]);
f = 2.*cos(y(:,1)) + y(:,2);
plot(t,f)

Кроме того, я хочу порекомендовать вам не использовать global. Вам это здесь не нужно, вы можете определить omega в качестве входного аргумента для modelo10: 

function dy = modelo10(t,y,omega)
   dy = zeros(3,1);
   dy(1) = y(2) - y(3);
   dy(2) = -3*y(2) - 5*sin(omega*y(1));
   dy(3) = y(1)*y(2);
end

и затем позвоните ode45 следующим образом: 

[t,y] = ode45(@(t,y)modelo10(t,y,omega), [0,20], [0,4,1]);

Здесь @(t,y)modelo10(t,y,omega) определяет анонимную функцию, которая несет в себе значение omega. Эта анонимная функция имеет два входных аргумента (t и y), как того требует ode45.

Наконец, вы можете упростить свой код, определив modelo10 в одной строке: 

function dy = modelo10(t,y,omega)
   dy = [ y(2) - y(3); -3*y(2) - 5*sin(omega*y(1)); y(1)*y(2) ];
end

Теперь вы можете сделать: 

modelo10 = @(t,y) [ y(2) - y(3); -3*y(2) - 5*sin(omega*y(1)); y(1)*y(2) ];
[t,y] = ode45(modelo10, [0,20], [0,4,1]);
2 голосов
/ 06 мая 2019 
global omega;
omega = 3;

[T,Y] = ode45(@modelo10, [0,20], [0,4,1]);

  • Отсюда Y является функцией T -> Нет необходимости определять новую функцию Y = Y (T) для T = [0,20], это больше не переменная, а действительные числа Y = [Y1 (0), ..., Y1 (20); Y2 (0), ..., Y2 (20); Y3 (0), ..., Y3 (20)]

  • Я назначаю функцию переменной f 

    f = @(t) 2.*cos(y(:,1)) + y(:,2) --> f = abs( 2.*cos(Y(:,1)) + Y(:,2));

    нижний регистр y и верхний регистр Y - две разные переменные. Функция абсолютного значения в Matlab - это abs

  • И, наконец, нанесите на график значения t 

    Use plot instead of fplot

    На самом деле вам не нужно указывать входное значение для функции f, поскольку оно имеет был уже включен при расчете ode45, тогда как fplot берет ключ аргумент, а именно, выражение функции 

    (function handle --> function defined as follow f = @(x)_____).

    Входные данные могут быть дополнительно добавлены в качестве второго аргумента.

    Это необязательно; по умолчанию fplot оценивает функцию в диапазоне [0 5]

  • Здесь f - действительное число, а не выражение функции, оно уже вычислено 

        plot(T, f);

enter image description here Функция f График для T в диапазоне от 0 до 20

...