RK4 Method for Windkessel error

I have these errors

Attempted to access P(1.16667); index must be a positive integer or logical.

Error in @(t,P,Q)(-P(t)/(R*C))+(Q(t)/C)

Error in RK4_WINDKESSEL (line 24)

k2=f(t(i)+0.5*k1*h, P(i)+0.5*k1, Q(i));

clear all

close all

clc

%parameters

R=1.5; %resistance

C=4; %compliance

h=1; %step size

t=1:h:100; %time vector

P(1)=1; %initial conditions

Q(1)=1; %initial conditions

f=@(t,P,Q) (-P(t)/(R*C))+(Q(t)/C);

%RK4 LOOP

for i=1:ceil(100/h)

if (i>=1 || i<=5)

Q(i)=i+1; %ramp

end

Q(i+5)=Q(i); %ramp

t(i+1)=t(i)+h;

k1=f(t(i),P(i),Q(i));

k2=f(t(i)+0.5*k1*h, P(i)+0.5*k1, Q(i));

k3=f(t(i)+0.5*k2*h, P(i)+0.5*k2);

k4=f(t(i)+k3*h, P(i)+0.5*k3);

P(i+1)=P(i)+(1/h)*6*(k1+2*k2+2*k3+k4);

end

plot(t,P)

xlabel('time')

ylabel('Pressure')

Answer 1

[James Tursa]

Multiple errors. When you call your function handle f, it is assumed that the P and Q are values passed at time t. They are not functions or function handles. So just use P and Q, not P(t) and Q(t):

f=@(t,P,Q) (-P/(R*C))+(Q/C);

Then in your RK4 code, the k's are not part of the t calculations. So the calls should look like this for the t part:

k1=f(t(i)      , etc);
k2=f(t(i)+0.5*h, etc);
k3=f(t(i)+0.5*h, etc);
k4=f(t(i)+    h, etc);

For the etc parts, I am not sure how to advise yet. Do you have two different state variables P and Q? And both have time derivatives that you know the equations for? If so, then you need to rewrite f and the k's code to account for this. Maybe separate f's and k's for the P and Q variables, or maybe carry a 2-element state vector instead.

Can you post the differential equation(s) that you are working with? Then we can advise further.

[oumi k]

Thank you for your guidance. The main differential equation is as follows:

The conditions for Q(t) as it is a ramp function which is an input in the differential equation.

[James Tursa]

This partially clears things up. You only have one state variable, P, and the Q is an input forcing function.

But the ramp/sawtooth function Q(t) looks strange. You have the basic ramp defined over a range of 0-4, but then repeats starting at 5. What happens between 4 and 5? Also the use of i in these equations is confusing, since it is being used both as a subscript and as a value of the function itself. Is this supposed to be a sawtooth function? This isn't well defined and needs to be explained.

[oumi k]

Yes, it is a sawtooth function. The i index was used for defining the sawtooth min and max range so a repetition from 1-5 on the vertical axis and the horizontal axis would be continuous time t.