Using embedded matlab funciton to find high accuracy derivative

Hi, I want to use embedded matlab function to take a high accuracy derivative of a sin wave. I have used sine wave block from simulink library, with 101 samples per period, 0 offset, 1 amplitude and sample time of 0.01. The solver is configured as fixed step solver (ode 4) with the same step size of 0.01. The simulation time is set as 1 sec.

The sine signal is named as B and here is the code I am using in the embedded funciton block. The formula for the derevative that I want to use uses central difference technique and is as follows.

B'(n))=(-B(n+2)+8B(n+1)-8B(n-1)+B(n-2))/(12.*step size); where, B is the input signal, n is the current time step I have successfuly tested the formula using command line and script and it is working, but when I use the same code in embedded matlab function, I get out of bounds error.

Here is the code that I am using:

function D = fcn(B) D=zeros (101,1); for n=3:99 D(n)=(-B(n+2)+8.*B(n+1)-8.*B(n-1)+B(n-2))./(12.*0.01); end This is the error message that I receive: Runtime error: Index into array out of range Attempted to access 5 element of data B. The valid index range is 1 to 1. I think it is due to the signal size problem, but I dont know how to fix it. Please help me, I will really appreciate you help! Thanks

Answer 1

[Titus Edelhofer]

Hi,

this is a misunderstanding of how Simulink works: Simulink computes every timestep on it's own. So when your embedded MATLAB function block get's called the input is the input from one single step (index range is 1, as the error says), not all 101 values.

So: for the time step at time T you don't have the input values of T+h, T+2h! These are values from the "future". You will have to use backward differencing (and use some default values for the history for T=0, T=h) if you want to use the embedded MATLAB function block.

Titus

[Amir Khan]

Hi Titus,

Thanks for your help!I got your point now. The formula for backward differencing in terms of D and B would be something like:

D(n)=3B(n)-4B(n-1)+B(n-2)/2(0.01)..

Assume that at T+2h the value is 1 and at T+2h, the value is 0. What would be the proper code for it then? I am fairly new to simulink and am struggling with it. Thanks a lot for you help!

[Titus Edelhofer]

Hi Amir,

you will use UnitDelay blocks to get the output of the sine at T-h, T-2h (just add another UnitDelay). The have an initial condition (which you probably best leave as zero).

Titus

[Amir Khan]

Thanks Titus! Really appreciate your help!