Question

Write MATLAB scripts to solve problem 5.21 from the textbook (bisection, specific heat) four times:

a) with the bisection method (starting interval [0, 1200])

b) with the false-position method (starting interval [0, 1200])

c) with the Newton-Raphson method (x0 = 0)

d) with the secant method ((x0, x1) = (0, 1))

use percent relative error = 10^-6 as stopping criterion.

5.21-Mechanical engineers, as well as most other engineers, use thermodynamics extensively in their work. The following polynomial can be used to relate the zero-pressure specific heat of dry air c_p kJ/(kg K) to temperature (K):

c_p = 0.99403 + 1.671 × 10^?4T + 9.7215 × 10^?8T² ?9.5838 × 10^?11T³ + 1.9520 × 10^?14T ⁴

Develop a plot of c_p versus a range of T = 0 to 1200 K, and then determine the temperature that corresponds to a specific heat of 1.1 kJ/(kg K).

Answer 1

Matlab code:txt format

1.bisection.m

%%----Part 1:Bisection Method--------%%
%to apply Bisection Method f(a)*f(b)<0 for interval [a b]
%%----------------Given input------------%%
%%given interval:[a b]= [0 1200]
clear all
clc
a=0;
b=1200;
%%tolerance iput means up to what % of accuracy in y
tor=1e-6;
%%-----Bisection Method iteration---------%%
%%checking weather f(a).f(b)<0 or >=0
if(myfun(a)*myfun(b)>=0)
error('choose another interval for roots finding');
else%if f(a).f(b)<0 then calculate c middle of a and b
c=(a+b)/2;
x(1)=c;
y(1)=myfun(c);
yval=abs(myfun(c));%calculating the value of y at x=c here myfun
%is a function that calulate the value of y
s=0;
while yval>tor
  
s=s+1;
  
if (myfun(a)*myfun(c)<0)
b=c;
else
a=c;
end
c=(a+b)/2;
x(s+1)=c;
y(s+1)=myfun(c);
yval=abs((y(s+1)-y(s)));
end
end
%%---------roots printing--------%%
fprintf('the approximate Temperature is:%.4f',c);

Output:

the approximate Temperature is:544.0842

2.regula.m

%%----Part 2:Regula Falsi Method--------%%
%to apply Regula Falsi Method f(a)*f(b)<0 for interval [a b]
%%----------------Given input------------%%
%%given interval:[a b]= [0 1200]
clear all
clc
a=0;
b=1200;
tor=1e-6;%solution error tolerance

%%-----Secant Method iteration---------%%
if(myfun(a)*myfun(b)>=0||(myfun(b)-myfun(a))==0)
error('choose another interval for roots finding');
else
c=(myfun(b)*a-myfun(a)*b)/(myfun(b)-myfun(a));
x(1)=c;
y(1)=myfun(c);
yval=abs(myfun(c));
ii=0;
while yval>tor
ii=ii+1;
if (myfun(a)*myfun(c)<0)
b=c;
c=(myfun(b)*a-myfun(a)*b)/(myfun(b)-myfun(a));
else
a=c;
c=(myfun(b)*a-myfun(a)*b)/(myfun(b)-myfun(a));
end
x(ii+1)=c;
y(ii+1)=myfun(c);
yval=abs(((y(ii+1)-y(ii))));
end
end
%%--------Print roots----%%
fprintf('the approximate Temperature is:%.4f\n',c);

Output:

the approximate Temperature is:544.0875
3.newton.m

%%----Part 3:Newton Raphson Method--------%%
%no such condition that f(a)*f(b) should be less than or greater zero

%%----------------Given input------------%%
%%given interval [0 1200]
a=0;
%enter the value of accuracy level
accur=1e-6;
  
%---------Newton Raphson Iteration---------%%
if(myfundiff(a)==0)
error('choose another intial guess for roots finding');
else
r=a-(myfun(a))/(myfundiff(a));
yval=abs(myfun(a));
ik=0;
while yval>accur
ik=ik+1;
a=r;
r=a-(myfun(a))/(myfundiff(a));
yval=abs(myfun(a));
end
end
%----------Printing roots------%
fprintf('the approximate Temperature is:%.4f',r);

Output:

the approximate Temperature is:544.0875

4.secant.m

%%----Part 3:Secant Method--------%%
%no such condition that f(a)*f(b) should be less than or greater zero

%%----------------Given input------------%%
%%given interval [0 1200]
clear all
clc
a=0;
b=1200;
tor=1e-6;%solution error tolerance;

%%-----Secant Method iteration---------%%
if((myfun(b)-myfun(a))==0)
error('choose another interval for roots finding');
else
c=(myfun(b)*a-myfun(a)*b)/(myfun(b)-myfun(a));
yval=abs(myfun(c));
ij=0;
while yval>tor
ij=ij+1;
a=b;
b=c;
c=(myfun(b)*a-myfun(a)*b)/(myfun(b)-myfun(a));
yval=abs(myfun(c));
end
end
%%--------Print roots----%%
fprintf('the approximate Temperature is:%.4f\n',c);

Output:

the approximate Temperature is:544.0876

Other two supporting m file for above Matlab code:

1.mufun.m

%here x and y stand for T and Cp respectively
function y=myfun(x)%here function name is myfun

y=-1.1+0.99403+(1.671e-4)*x+(9.7215e-8)*x^2-(9.5838e-11)*x^3+(1.9520e-14)*x^4;

end

2.myfundiff.m

%here x stands for T
function diff_f=myfundiff(x)%here function name is myfundiff

diff_f=(1.671e-4)+2*(9.7215e-8)*x-3*(9.5838e-11)*x^2+4*(1.9520e-14)*x^3;

end

Note: Put all the above Matlab file in same folder then run

Matlab code for Plot of Cp vs T:

plotcp_temp.m

%%------Graph Cp vs T-------%%
x=0:10:1200;
y=0.99403+(1.671e-4).*x+(9.7215e-8).*x.^2-(9.5838e-11).*x.^3+(1.9520e-14).*x.^4;
plot(x,y)
xlabel('T(in K)')
ylabel('Cp(in Kj/kgK)')

Output:

the value of temperature T for which Cp=1.1 Kj/kgK is T=545K (From Graph)