Same time steps for ODE function and OutputFcn

조회 수: 5 (최근 30일)
Amavi Silva
Amavi Silva 2021년 8월 19일
답변: Alan Stevens 2021년 8월 19일
In the following model, I need to know the values of the model parameter 'k' along with 'y'. To find k, I used 'OutputFcn'. However, after runing the model, I get different lengths for the k array (13) and 'y' array (94). But what I need to know is the value of k each iteration uses to compute each 'y' value (or simply, similar lengths for both arrays). In addition, I see that 'x' (equivalent to time) called in the ODE function also has a different array length (94) compared to that of the same 'x' called in OutputFcn (13). I get that the reason behind this is the fact that OutputFcn generates parameter values after each integration timestep rather than after every call of the ODE function. I thought I will perhaps solve the issue by asssigning fixed time steps to the ODE function. But the issue remains the same. Is there anyway I can modify the script where both OutputFcn and ODE function would return the same array lengths for 'k' and 'y' using same 'x'? Or else is there any other method I can get my expected result? Any help is much appreciated.
Thank you!
global k_all kx
options = odeset('Reltol', 1e-3, 'OutputFcn', @get_process, 'MaxStep', 100, 'InitialStep', 1);
kx = []; % both l and x values together
k_all = []; % All values
t0 = 0;
tf = 2000;
tspan = t0:0.005:tf;
[x,y] = ode45(@myode,tspan,1,options);
function dydx = myode(x, y)
global kx
k = x.^2 + y.^2;
kx = [k x];
dydx = x + k.*y;
end
function [status] = get_process(x, y, flag);
global kx k_all
k_all = [k_all; kx];
status = 0;
end

답변 (1개)

Alan Stevens
Alan Stevens 2021년 8월 19일
You could simply calculate
k = x.^2 + y.^2;
immediately after
[x,y] = ode45(@myode,tspan,1,options);
However, have you really looked at your ode? You have, in effect,
dy/dx = x + x.^2.*y + y.^3;
and you start from y = 1.
y blows up extremely rapidly and the ode solvers won't converge after x is slightly greater than 0.4, so there's no way you are going to get to x = 2000.

카테고리

Help CenterFile Exchange에서 Ordinary Differential Equations에 대해 자세히 알아보기

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by