Problems with a Newton-Raphson code
1 view (last 30 days)
Show older comments
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
0 Comments
Answers (2)
Brendan Hamm
on 6 Apr 2015
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.
0 Comments
Jan
on 6 Apr 2015
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.
0 Comments
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!