Step response of nonlinear differential equation

Hello!

I've got a following non-linear differential equation: y_dot = (A1/B)*(x^2)*y + (A2/B)*x*y + (A3/B)*y + C

• y_dot = (A1/B)*(x^2)*y + (A2/B)*x*y + (A3/B)*y + C

Where y is an output I'm interested in and x is input. Is it possible to get a step response for such equation? I've already got it implemented in Simulink, but scopes don't show things like settling time and other statistics I could use to derive PID controller parameters.

Answer 1

[Vijay Periasamy]

I have written a function for this. Step_ODE(fhan, Solver, t_s, t_t, Val_ini, Val_fin, ini) Please have a look at: https://au.mathworks.com/matlabcentral/fileexchange/65201-step-ode-fhan--solver--t-s--t-t--val-ini--val-fin--ini-

[Vijay Periasamy]

Please find updated link below: https://au.mathworks.com/matlabcentral/fileexchange/65201-step-response-of-a-non-linear-differential-equation-system

Answer 2

[Vijay Periasamy]

If you could provide details on why it is not working, I can help. I am pasting the example below to answer some of your questions.

------------------------------ Description of input arguments: ------------------------------ fhan - Function handle for the differential equation function Solver - String name of the ODE solver t_s - Step time t_t - Total simulation time Val_ini - Vector of initial values (before step) for the inputs or parameters Val_fin - Vector of final values (after step) for the inputs or parameters ini - Intial value vector for the differential equations ------------------------------ Example usage for a 2 state, 2 parameter system: Val_ini = [0.3 0.4]; Val_fin = [0.5 0.1]; t_s = 150; t_t = 300; ini = [0 0]; Solver = 'ode23t'; fhan = @myfun; [t,y] = Step_ODE(fhan, Solver, t_s, t_t, Val_ini, Val_fin, ini); plot(t,y(:,1)) figure; plot(t,y(:,2))

function dydt = myfun(~,y, Lam)

dydt = zeros(2,1); dydt(1) = -Lam(1) * y(1) + 5 ; dydt(2) = -Lam(2) * y(2) + 10 * y(1);

end ------ Firstly, the function needs to be in the working directory. A function like myfun above needs to be created where you input your differential equations. How many inputs or parameters that needs stepping do you have in your model? The above example has 2: Lam(1) and Lam(2). In above example Lam(1) is changed from 0.3 to 0.5 and Lam(2) is changed from 0.4 to 0.1 at 150 seconds. See below. Val_ini = [0.3 0.4]; Val_fin = [0.5 0.1]; This is how you change the step magnitude (initially 0.3 to 0.5; initially 0.4 to 0.1) It is not necessary that both have to be changed at same time. You can just change one. For example; Val_ini = [0.3 0.4]; Val_fin = [0.3 0.1]; Here only Lam(2)is changed. What do you mean by range of step? If you mean step time, then it can be done as follows. For example you need the step to be at 100 seconds instead of 150 seconds, then instead of t_s = 150; use t_s = 100;

Initial and final vectors consist of the initial and final values of the input you want to step. They are vectors because there can be more than one input you are interested in. In case, you have only one input, they need not be vectors and can be like the following. Val_ini = 0.3; Val_fin = 0.5;

[sangram dogra]

okay thanks for reply but the error it shows is that too many input arguments. I am giving the function code. in this my input is f(1) and f(2). f(2), i am keeping constant, and the initial value for my function are ini = [5.37545 2.2433*10^-1 3.1308*10^-3 6.2616*10^-1].

if true
  % code
%%where po here is
po=(2*f0*y(2)*zl*exp(-el/(r*t))/(zd*exp(-ed/(r*t))+zc*exp(-ec/(r*t))))^0.5;
%%here goes the code
%y(1)=Cm; y(2)=C1; y(3)=D.; y(4)=D1;
dy=zeros(4,1);
dy(1)=-(zp*exp(-ep/(r*t))+zm*exp(-em/(r*t)))*y(1)*po-f(2)*y(1)/v+f(2)*cm/v;
dy(2)=-zl*exp(-el/(r*t))*y(2)-(f(2)*y(2)/v)+f(1)*cl/v;
dy(3)=(0.5*zc*exp(-ec/(r*t))+zd*exp(-ed/(r*t)))*po^2+zm*exp(-em/(r*t))*y(1)*po-f(2)*y(3)/v;
dy(4)=m*(zp*exp(-ep/(r*t))+zm*exp(-em/(r*t)))*y(1)*po-f(2)*y(4)/v;
end

[Vijay Periasamy]

That error can occur if you passed a function more input arguments than it expected to receive. Since your function has 2 inputs and 4 variables, your call to the function should look like the following. Val_ini should be a vector of length 2, with your desired intial values of f(1) and f(2). Val_fin should be a vector of length 2, with your desired final values of f(1) and f(2).

Val_ini = [0.3 0.4]; 
Val_fin = [0.5 0.1]; 
t_s = 150; 
t_t = 300; 
ini = [5.37545 2.2433*10^-1 3.1308*10^-3 6.2616*10^-1]; 
Solver = 'ode23t'; 
fhan = @myfun; 
[t,y] = Step_ODE(fhan, Solver, t_s, t_t, Val_ini, Val_fin, ini);

And the function should look like this. f

unction dy = myfun(~,y, f) 
po=(2*f0*y(2)*zl*exp(-el/(r*t))/(zd*exp(-ed/(r*t))+zc*exp(-ec/(r*t))))^0.5;
%%here goes the code
%y(1)=Cm; y(2)=C1; y(3)=D.; y(4)=D1;
dy=zeros(4,1);
dy(1)=-(zp*exp(-ep/(r*t))+zm*exp(-em/(r*t)))*y(1)*po-f(2)*y(1)/v+f(2)*cm/v;
dy(2)=-zl*exp(-el/(r*t))*y(2)-(f(2)*y(2)/v)+f(1)*cl/v;
dy(3)=(0.5*zc*exp(-ec/(r*t))+zd*exp(-ed/(r*t)))*po^2+zm*exp(-em/(r*t))*y(1)*po-f(2)*y(3)/v;
dy(4)=m*(zp*exp(-ep/(r*t))+zm*exp(-em/(r*t)))*y(1)*po-f(2)*y(4)/v;
end

[Vijay Periasamy]

you mentioned f(2) is constant. Then why dont you declare it as a constant like variables f0, zp, zc etc? Perhaps you want to step it later. In the following example I keep f(2) constant with value 1 and f(1) stepping from value 1 to 3.

Val_ini = [1 1]; 
Val_fin = [1 3];

[Vijay Periasamy]

You did not mention the initial and final values of f(1). I assume them to be 1 and 1.5 and got following results.

Val_ini = [1 1]; 
Val_fin = [1.5 1]; 
t_s = 150; 
t_t = 300; 
ini = [5.37545 2.2433*10^-1 3.1308*10^-3 6.2616*10^-1]; 
Solver = 'ode23t'; 
fhan = @myfun; 
[t,y] = Step_ODE(fhan, Solver, t_s, t_t, Val_ini, Val_fin, ini); 
figure
plot(t,y(:,1)) 
figure; 
plot(t,y(:,2)) 
figure; 
plot(t,y(:,3)) 
figure; 
plot(t,y(:,4)) 

function dy = myfun(~,y,f) 
cl=8;
cm=6;
ec=2.9442*10^3;
ed=2.9442*10^3;
em=7.4478*10^4;
el=1.2550*10^5;
ep=1.8283*10^4;
f0=.58;
m=100.12;
r=8.314;
t=335;
zc=3.8223*10^10;
zd=3.1457*10^11;
zm=1.0067*10^15;
zl=3.7920*10^18;
zp=1.7700*10^9;
v=0.1;

po=(2*f0*y(2)*zl*exp(-el/(r*t))/(zd*exp(-ed/(r*t))+zc*exp(-ec/(r*t))))^0.5;
dy=zeros(4,1);
dy(1)=-(zp*exp(-ep/(r*t))+zm*exp(-em/(r*t)))*y(1)*po-f(2)*y(1)/v+f(2)*cm/v;
dy(2)=-zl*exp(-el/(r*t))*y(2)-(f(2)*y(2)/v)+f(1)*cl/v;
dy(3)=(0.5*zc*exp(-ec/(r*t))+zd*exp(-ed/(r*t)))*po^2+zm*exp(-em/(r*t))*y(1)*po-f(2)*y(3)/v;
dy(4)=m*(zp*exp(-ep/(r*t))+zm*exp(-em/(r*t)))*y(1)*po-f(2)*y(4)/v;
end
and got following results.
<</matlabcentral/answers/uploaded_files/129235/y1.bmp>>

<<

>>

Answer 3

[Sherin Ashraf-Hanna]

𝑓(𝑡) =

𝑓(0)𝑒^(k*t)

[Sherin Ashraf-Hanna]

applied to 5730 years, 77.45%intials, t=0, 14 is present how to apply this on matlab throught the formula I have above?

thank you