Problems with a Newton-Raphson code

function pnew = newtraph(p4,p1,d,p,e,f)

p1=20;

p4=500;

k1=1.4;

k4=1.4;

r1=287;

r4=287;

c1=348.92;

c4=348.92;

p=1;

dp=.0000001;

d=((k4-1)/(2*k1))*(c1/c4)

e= (k1+1)/(2*k1)

f= (2*k4)/(k4-1)

while abs(f(p)) > 1e-6

fp =(p4/p1)*(1-((d*(p-1)/(1+e*(p-1))))^f)-p

fpmdp= (p4/p1)*(1-((d*((p-dp)-1)/(1+e*((p-dp)-1)))))^f-(p-dp)

fppdp= (p4/p1)*(1-((d*((p+dp)-1)/(1+e*((p+dp)-1)))))^f-(p+dp)

dfdp= (fppdp-fpmdp)/2*dp

pnew=fp-(fp/dfdp)

p=pnew

end

I am trying to use a newton raphson method to find pnew for a shock tube in my gas dynamics class. I get the error message below.

Subscript indices must either be real positive integers or logicals.

Error in Untitled (line 34) end

With the given values from above, the iteration should produce a pnew=4.0471

Any help would be much appreciated, thanks

Answer 1

[Brendan Hamm]

You have a line:

while abs(f(p)) > 1e-6

which is an indexing into the variable f , with the indexing variable p. But p changes on every iteration in the while loop. It changes to the value 9.230769236418623e+13 on the first loop and you get an error when you try and use this value as an index into f at the while loop line.

Answer 2

[Jan]

The debugger helps you to indentify the source of problems. Set a breakpoint in your code and step through the function line by line. Then you will find out, what Brendan has mentioned already.