# Computational Physics by 2lROJe5M

VIEWS: 0 PAGES: 33

• pg 1
```									Computational Physics
Chapter 3

Problem 3.23, 3.25, 3.30
Solution
The solution is the root of the function f(T)=100(1+AT-BT2)-200
A plot of the function for the domain[100,500] is obtained with MATLAB by typing:

fplot(‘100*(1+3.90802E-3*x-0.580195E-6*x^2)-200’,[100,500])
(a) Use the user-defined function BisectionRoot(in 3.14) with 0.001 for TolMax

f(T)=100(1+AT-BT2)-200

function f=FunHW3_23(T)
A=3.90802E-3;B=0.580195E-6;
f=100*(1+A*T-B*T^2)-200;

In the command window, use the BisectionRoot to obtain the root

>>BisectionRoot('FunHW3_23', 250,300, 0.001)

Function of which          Domain Tolerance
we want to obtain the root
(b) Use MATLAB’s built-in fzero function

In the command window,

>>fzero('FunHW3_23',250)

A value of x near to where
The function to be solved
the function crosses the axis
Solution
The solution is the root of the function f(Ts)=πdL[h(Ts-Tair)+εσSB(T4s-T4surr)]-18405
A plot of the function for the domain[0,1000] is obtained with MATLAB by typing:

fplot('pi*0.1*25*(10*(T-298)+0.8*5.67E-8*(T^4-298^4))-18405',[0,1000])
(a) Use the user-defined function BisectionRoot. that was developed in Problem
3.14 with a starting interval of [0,100]. For maximum tolerance use 10-2.

f(Ts)=πdL[h(Ts-Tair)+εσSB(T4s-T4surr)]-18405

function f=FunHW3_25(T)
L=25; d=0.1; e=0.8; s=5.67E-8; Q=18405;
h=10; Ta=298;
f=pi*d*L*(h(T-Ta)+e*s*(T^4-Ta^4))-Q;

In the command window, use the BisectionRoot to obtain the root

>>Xs=BisectionRoot('FunHW3_25',0,1000,0.01)

(b) Using MATLAB's built-in fzero function(in the Command Window) to find the
root:
In the command window,

>>fzero('FunHW3_25',400)
Solution

function y=Fun3_30(x)
q=1.6022E-19; k=1.3806E-23; Voc=0.5;
T=297;
qkT=q/(k*T);
y=exp(qkT*x)*(1+qkT*x)-exp(qkT*Voc);
The function is used in the following program(script file) to make the plot

V=0.1:0.01:0.45;
n=length(V);
for i=1:n
f(i)=Fun3_30(V(i));
end
plot(V,f)
xlabel('Vmp'); ylabel('f')

The figure shows the root is between 0.4 and 0.45
(a) To use the fixed-point iteration method the equation e^(...)=e(...) is solved for
Vmp:

This iteration function is used in the following MATLAB program. For starting
point the value Vmp=0.5V is used.
The iterations are terminated by using Eq(3.18) with epsilon=0.001

%Fixed Point iteration method
clear all
q=1.6022E-19; k=1.3806E-23; Voc=0.5; T=297;
qkT=q/(k*T);
x=0.5
for i=1:7
xnext=(1/qkT)*log(exp(qkT*Voc)/(1+qkT*x))
if abs((xnext-x)/x)<0.001
break
else
x=xnext;
end
end
(b) MATLAB's built-in function fzero is used in the Command window to find the
root

>>Xs=fzero('Fun3_30',0.5)
Chapter 4

Problem 4.33, 4.34, 4.41
Solution

Using MATLAB’s left division command
>> A=[0.2 0.2 0.3 0.2 0.2; 28 1 0 0 0.1; 0 18 12 2.4 16; 0 0 10 0 1; 0 0 0 10 2; 0 0
0 0 18];
>>b=[3.4;20.5;170;49;39.8;96.3];
>>x=A\b

Thus, the concentrations of the different species in the sample are:
.
nCH4   = 0.6634 nC2H4 = 1.3909 nC2H6 = 4.3650
nC H
3   6   = 2.9100 nC3H8 = 5.3500
function x = GaussJordan(a,b)
% The function solve a system of linear equations ax=b using the Gauss
% elimination method with pivoting. In each step the rows are switched
% such that pivot element has the largest absolute numerical value.
% Input variables:
% a The matrix of coefficients.
% b A column vector of constants.
% Output variable:
% x A column vector with the solution.

ab = [a,b];
[R, C] = size(ab);
for j = 1:R
% Pivoting section starts
pvtemp=ab(j,j);
kpvt=j;
% Looking for the row with the largest pivot element.
for k=j+1:R
if ab(k,j)~=0 & abs(ab(k,j)) > abs(pvtemp)
pvtemp=ab(k,j);
kpvt=k;
end
end
% If a row with a larger pivot element exists, switch the rows.
if kpvt~=j
abTemp=ab(j,:);
ab(j,:)=ab(kpvt,:);
ab(kpvt,:)=abTemp;
end
% Pivoting section ends
ab(j,:)= ab(j,:)/ab(j,j);
for i = 1:R
if i~=j
ab(i,j:C) = ab(i,j:C)-ab(i,j)*ab(j,j:C);
end
end
end
x=ab(:,C);
(a) Solve the system of equations using the user-defined function GaussJordan
developed in Problem 4.24.

% Solution of HW4_34a, script file
clear all
A=[-1 0 -0.7071 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 -1 -0.7071 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0.7071 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 -1 0 -0.7071 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 -1 0 0 0.7071 1 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 -1 -0.7071 0 0 1 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0.7071 0 1 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 -1 0 -0.7071 0 0 0.9806 0.7071 0 0 0 0 0 0
0 0 0 0 0 0 0 -1 0 0 0.7071 0 1 0.1961 0.7071 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 -1 -0.7071 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0.7071 1 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 1 0 0 0 0 0
000000000000100000000
0 0 0 0 0 0 0 0 0 0 0 0 0 -0.9806 0 0 0 0.9806 0.7433 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 -0.1961 0 0 1 0.1961 0.669 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.7071 -1 0 0 0 1 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.7071 0 -1 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.7433 -1 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.669 0 -1
0 0 0 0 0 0 0 0 0 0 0.9806 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.1961 0 0 1];
b=[0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 8000; 0; 0; -5000; 0;
1000*cosd(60); 10000*sind(60)];
c = GaussJordan(A,b)
(b) Solve the system of equations using MATLAB’s left division operation.

% Solution of HW4_34b, script file
clear all
A=[-1 0 -0.7071 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 -1 -0.7071 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0.7071 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 -1 0 -0.7071 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 -1 0 0 0.7071 1 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 -1 -0.7071 0 0 1 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0.7071 0 1 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 -1 0 -0.7071 0 0 0.9806 0.7071 0 0 0 0 0 0
0 0 0 0 0 0 0 -1 0 0 0.7071 0 1 0.1961 0.7071 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 -1 -0.7071 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0.7071 1 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 1 0 0 0 0 0
000000000000100000000
0 0 0 0 0 0 0 0 0 0 0 0 0 -0.9806 0 0 0 0.9806 0.7433 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 -0.1961 0 0 1 0.1961 0.669 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.7071 -1 0 0 0 1 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.7071 0 -1 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.7433 -1 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.669 0 -1
0 0 0 0 0 0 0 0 0 0 0.9806 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.1961 0 0 1];
b=[0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 8000; 0; 0; -5000; 0;
1000*cosd(60); 10000*sind(60)];
c=A\b
(a) Determine the eigenvalues (frequencies), and the corresponding wavelengths
( λ=2πc/ω(where c=3 *108m/s is the speed of light).)
% Solution of Problem 4.41. Script file.
clear all
kCH=5.92e2; kCC=15.8e2; mH=1.6605e-27; mC=12*1.6605e-27;
c=3e8; M=zeros(4);
M(1,1)=kCH/mH; M(1,2)=-M(1,1);
M(2,1)=-kCH/mC; M(2,2)=(kCH+kCC)/mC; M(2,3)=-kCC/mC;
M(3,2)=M(2,3); M(3,3)=M(2,2); M(3,4)=M(2,1);
M(4,3)=M(1,2); M(4,4)=-M(4,3);
[evectors values]=eig(M); eigenvalues=diag(values);
temp=sqrt(eigenvalues); icount=1;
% Looking for the real eigenvalues and corresponding eigenvectors.
for i=1:length(temp)
if imag(temp(i))==0
number=temp(i);
omega(icount)=number;
eigenvectors(:,icount)=evectors(:,i);
icount=icount+1;
end
end
disp('The real eigenvalues are:')
omega
disp('The corresponding wavelengths (except for the zero eigenvalue) in meters are:')
for i=2:length(omega)
lambda=2*pi*c./omega;
end
lambda
disp('The eigenvectors corresponding to the eigenvalues are:')
eigenvectors
(b) Determine the eigenvectors corresponding to the eigenvalues found in part (a).
From the eigenvectors, deduce the relative motion of the atoms (i.e., are they
moving towards, or away from each other?)

The columns of the eigenvector matrix are the individual eigenvectors
themselves. They represent amplitudes of the motion or vibration of each atom
in the C2H2 molecule.
The first column represents the eigenvector for the eigenvalue associated with
the symmetric stretch of the molecule where the hydrogen atoms at the end and
their neighboring carbon atoms are moving apart. This motion corresponds to a
wavelength of λ=5.05mm in the infrared portion of the electromagnetic spectrum.
The second column indicates the leftmost H atom moving to the right, the
leftmost C atom moving to the left, the rightmost C atom moving to the right and
the rightmost H atom moving to the left. In other words, the H-C and C-H bonds
are vibrating such that they are being compressed and the C atoms in the interior
are being stretched apart. This corresponds to a wavelength of λ=2.96mm , in
the near-infrared.
The third column is the eigenvector corresponding to the eigenvalue that
represents an asymmetric stretch of the molecule. The two C atoms are moving
to the left while both H atoms are moving to the right. Thus, the left side of the
molecule is being compressed while the right side of the molecule is being
stretched. This corresponds to a wavelength of λ=3.03mm , also in the infrared.
.
Chapter 5

Problem 5.24, 5.27
Solution
To solve the problem the equation k=bem(d/D)is written in the form:
ln(k)=m(d/D)+ln(b)
Then, linear least-squares regression is used for finding the coefficients m and
K that best fit the data. (See equations in the second row of Table 5-2.)
The following program written in a script file determines the constants and plot
the data points and the curve that best fit the data.

% Solution of HW5_24 script file
clear all
dD=[0.05 0.25 0.45 0.65 0.85];
k=[2.91 2.4 2.17 2.11 2.03];
n=length(k);
Y=log(k);
p=polyfit(dD,Y,1);
m=p(1)
b=exp(p(2))
x=linspace(dD(1),dD(n),30);
y=b*exp(m*x);
plot(x,y)
hold on
plot(dD,k,'*r')
xlabel('d/D'), ylabel('k')
legend('Best fit curve','Data points')
hold off
% Model prediction for d/D = 0.15
kFor015=b*exp(m*0.15)
m=
-0.4245
b=
2.7889
kFor015 =
2.6168
User defined function CubicSplines form Problem 5.22

function Yint=CubicSplines(x,y,xint)
% CubicSplines fits a set of cubic splines to n data points and then
% returns the interpolated value Yint at the desired coordinate Xint.
% Input variables:
% x A vector with the coordinates x of the data points.
% y A vector with the coordinates y of the data points.
% xint The x-value where interpolation is desired
% Output variables:
% Yint The interpolated y-value at x=xint
n=length(x); interval=1;
if n ~= length(y)
disp('ERROR: x and y do not have the same number of points');
stop
end
% calculate h_i
for i = 1:n-1
h(i) = x(i+1)-x(i);
end
%Start Thomas Algorithm
for i=2:n-2
b(i)=h(i);
end
b(1)=0;
for i=1:n-3
a(i)=h(i+1);
end
a(n-2)=0;
for i=1:n-2
d(i)=2*(h(i)+h(i+1)); r(i)=6*(((y(i+2)-y(i+1))/h(i+1))-((y(i+1)-y(i))/h(i)));
end
A=zeros(n-2,n-2);
for i=2:n-3
A(i,i)=d(i); A(i,i+1)=a(i); A(i,i-1)=b(i);
end
A(1,2)=a(1); A(1,1)=d(1); A(n-2,n-3)=b(n-2); A(n-2,n-2)=d(n-2);
a(1)=a(1)/d(1);
r(1)=r(1)/d(1);
for i=2:n-3
denom=d(i)-b(i)*a(i-1);
if(denom==0), error('zero in denominator'), end
a(i)=a(i)/denom;
r(i)=(r(i)-b(i)*r(i-1))/denom;
end
r(n-2)=(r(n-2)-b(n-2)*r(n-3))/(d(n-2)-b(n-2)*a(n-3));
ans(n-2)=r(n-2);
for i=n-3:-1:1
ans(i) = r(i) - a(i)*ans(i+1);
end
acoeff(1)=0; acoeff(n)=0;
for i=2:n-1
acoeff(i)=ans(i-1);
end
for j=1:n-1
if xint >= x(j) & xint <= x(j+1)
interval=j;
end
end
i=interval;
YintA=(acoeff(i)*((x(i+1)-xint)^3)/6/h(i));
YintB=(acoeff(i+1)*((xint-x(i))^3)/6/h(i));
YintC=((y(i)/h(i))-(acoeff(i)*h(i)/6))*(x(i+1)-xint);
YintD=((y(i+1)/h(i))-(acoeff(i+1)*h(i)/6))*(xint-x(i));
Yint=YintA+YintB+YintC+YintD;
(a) The following program in a script file solves the problem using the user-
defined function CubicSplines from Problem 5.22.

% Solution of HW5_27a script file
clear all
T=[5:2.5:30]*1000;
h=[3.3 7.5 41.8 51.8 61 101.1 132.9 145.5 171.4 225.8 260.9];
Tint=5000:500:30000;
n=length(Tint);
for i=1:n
hint(i)=CubicSplines(T,h,Tint(i));
end
plot(Tint,hint)
hold on
plot(T,h,'*r')
xlabel('Temperature (K)')
ylabel('Enthalpy (MJ/kg)')
legend('Interpolation','Data points',2)
text(20000,20,'Problem 5.27 Part (a)')
hold off
(b) The following program in a script file solves the problem using MATLAB’s
built-in function interp1l with the spline option.

% Solution of HW5_27b script file
clear all
T=[5:2.5:30]*1000;
h=[3.3 7.5 41.8 51.8 61 101.1 132.9 145.5 171.4 225.8 260.9];
Tint=5000:500:30000;
n=length(Tint);
hint=interp1(T,h,Tint,'spline');
plot(Tint,hint)
hold on
plot(T,h,'*r')
xlabel('Temperature (K)')
ylabel('Enthalpy (MJ/kg)')
legend('Interpolation','Data points',2)
hold off

```
To top