How can I record the values of pa and p0 in a table in this code?

while abs(y)<0.05
    
    d=(wn0^2/w^2)*(1-((l^2/pi^2)*((p035*Alum-N)/EI)));
    
    e=(wn0^2*l^2*pa35*Alum)/(w^2*pi^2*EI);
    
    [t,y]=ode45(@dydt,[0:1:2000],[0.01 0],[],d,e);
    
    y=y(length(t),1)
    
    if abs(y)<0.05
        p035=p035+h*133.332
        figure(1)
        hold on
        
    else
        plot(pa35/133.322,p035/133.322,'ro')
        pa35=pa35+h*133.332
        p035=0;
        y=0;
        if pa35>pacr35
            break
        end
    end

[M G]

Just to be clear, I need the values that are being ploted in this line: plot(pa35/133.322,p035/133.322,'ro')

Thank you!

Answer 1

[William Rose]

@M G, see attached. It runs. It produces the console output below. It also produces the plot below. There are comments in the code which explain what it is doing. Good luck with your work.

>> arteryOscillationWR

Enter Pressure step size in mmHg: 5

Calculating critical mean pressure: p0cr15=8000 Pa = 60.00 mmHg

Calculating critical pulse pressure: pacr15=8000 Pa = 60.00 mmHg

Calculating and plotting (Pa,P0) pairs

Pairs saved: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15.

Linear regression: intercept=58.6 mmHg, slope = -0.988.

>>

Please accept this answer.

[M G]

Thank you so much William!

Answer 2

[William Rose]

@Mustapha Gharib,

plot(pa35/133.322,p035/133.322,'ro')

The command above wil work if pa35 and p035 are vectors of the same length. If they are both length 1, then a single point will be plotted, with x=pa35 and y=p035. Division by 133.32 converts pressures from Pascals to mmHg.

Where do p035 and pa35 come from? They are not set in your code.

We cannot run your code, since you have not provided dydt().

You have a while loop with no end statement. Or the if-else has no end statement. Either way, you need another end statement.

Here is a section of code that needs attention. I have added comments.

%Next line returns a 2001x1 vector t, and 2001x2 array y. 
%We infer that y has 2 columns, because the initial conditions vector is 2x1.
%Does this ode45() command work?  The arguments at the end do not look normal
%or correct to me.
[t,y]=ode45(@dydt,[0:1:2000],[0.01 0],[],d,e); 
y=y(length(t),1)     %Keep first column of y, discard column 2.
%Does the next line do what you want?  y is a vector. 
%"if abs(y)<0.05" evaluates to True iff every element of y is <.05.
if abs(y)<0.05
    %Next line adds an offset to p035. h is in mmHg, p035 is in Pascals?
    p035=p035+h*133.332
    figure(1)
    hold on
else  %this option is chosen if any element of y>=.05
    %Next line plots pa35 versus p035.  Pressure versus pressure.    
    %Is that what you want? Where do vectors p035, pa35 come from?
    %They are not from the integration by ode45().
    plot(pa35/133.322,p035/133.322,'ro')
    %Next line adds an offset to pa35. h is in mmHg, pa35 is in Pascals?
    pa35=pa35+h*133.332
    p035=0;
    %Next line replaces the entire vector y with a single number: zero
    y=0;
    if pa35>pacr35
        break
    end
end

I hope that my comments above will help you evaluate and, if necessary, correct your code.