Numeric Methods to Solve ODEs(Ordinary Differential Equations )

Thank you very much Hamid. That is beautiful.

I used a different method - numerical integration (Tustin method) to do the same thing, which I think more accurate. Please note that the below code has not been tested.

FUNCTION_BLOCK FOP

TITLE = 'First Order Process'
// Y            A
//--- = --------------------
// U          Bs + C
//
//
// Use Tustin method to approximate
//
//      2     z-1
// s = ---*----------
//      T     z+1
//
//  Y               ATz + AT
// --- = ----------------------------------
//  U          (2B+CT)z + (CT-2B)
//
//
//
// (2B+CT)y(k+1) = ATu(k+1) + ATu(k) + (2B-CT)y(k)
//

VERSION: '0.0'
AUTHOR: ZLiang
NAME: FOP
FAMILY: PROCESS

// Block Parameters
VAR_INPUT
I_U:REAL;    //Input
I_A:REAL;   //System gain
I_B:REAL;   //Time constant
I_C:REAL:=1.0;  //Typical 1
I_Ts:TIME:=T#100MS;    //Sample time
END_VAR
VAR_IN_OUT
IQ_RST:BOOL:=FALSE; //Reset
END_VAR
VAR_OUTPUT
Q_Y:REAL; //Output
Q_ERR:BOOL; //Error
END_VAR
VAR
Y0:REAL:=0.0;   //Assume the initial condition is 0
Y1:REAL:=0.0;
U0:REAL:=0.0;   //Assume the initial condition is 0
Ts:REAL; //Sample time
END_VAR

// Instruction Section

Ts := (DINT_TO_REAL(TIME_TO_DINT (I_Ts)))/1000.0;

IF IQ_RST THEN
Y0:=0; //Assume the initial condition is 0
Y1:=0;
U0:=0.0;
Q_ERR:=FALSE;
IQ_RST:=FALSE;
ELSE
IF I_B < 0.005 AND I_B > -0.005 THEN
Q_ERR:=TRUE;
Y0:=0;
Y0:=0;
ELSE
// (2B+CT)y(k+1) = ATu(k+1) + ATu(k) + (2B-CT)y(k)
Y1 := ( ( I_A * Ts * I_U ) + (I_A * Ts * U0) + ( 2 * I_B - I_C * Ts ) * Y0 )/( 2 * I_B + I_C * Ts )
U0 := I_U;
Y0 := Y1;
Q_ERR := FALSE;
END_IF;
END_IF;

Q_Y := Y1;

END_FUNCTION_BLOCK

Cheers,

ZLiang

New question published by johnnybtb is split to a separate thread with the subject Numeric Methods to Solve ODEs(Ordinary Differential Equations ) in TIA.

Best regards
Jen_Moderator