RK4 Method for Windkessel error

1 vue (au cours des 30 derniers jours)
oumi k
oumi k le 21 Oct 2022
Commenté : oumi k le 22 Oct 2022
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')

Réponse acceptée

James Tursa
James Tursa le 21 Oct 2022
Modifié(e) : James Tursa le 21 Oct 2022
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.
  3 commentaires
James Tursa
James Tursa le 22 Oct 2022
Modifié(e) : James Tursa le 22 Oct 2022
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
oumi k le 22 Oct 2022
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.

Connectez-vous pour commenter.

Plus de réponses (0)

Tags

Community Treasure Hunt

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

Start Hunting!

Translated by