Figure 4.16:

Typical strain versus time data from a molecular dynamics simulation from data file \texttt {rohit.dat} on the website \jbrawweb .

Code for Figure 4.16

Text of the GNU GPL.

main.m


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
% jbr, 5/8/2007
% CBE 731, final exam problem

ndata = 150;
nplot = 200;
tplot = linspace(0, 2*pi, nplot)';
tdata = linspace(0, 2*pi, ndata)';

g = [3; 4];

% noise  is N(0,V) and V=vmax + (vmin-vmax)* (t/pi-1)^2
vmax = 3;
%vmin = 0.1;
%V=vmax + (vmin-vmax)*(tdata/pi-1).^2;
a = 5;
V=vmax*exp(-a*(tdata/pi-1).^2);

randn("seed", -1);

noise = sqrt(V).*randn(ndata,1);

X = [sin(tdata), cos(tdata)];
Xplot = [sin(tplot), cos(tplot)];
y = Xplot*g;
ymeas = X*g + noise;

gest = inv(X'*X)*X'*ymeas;
yest = Xplot*gest;
resid = ymeas - X*gest;
p = length(g);
samplevar = resid'*resid/(ndata-p);
alpha = 0.95;
Fstat = finv(alpha, p, ndata-p);
bbox = sqrt(2*p*Fstat*samplevar*diag(inv(X'*X)));

[gest, bbox]
plot (tplot, [y, yest], tdata, ymeas, 'o')
table = [tdata, ymeas];

% reduce precision in output "data"
[fid, msg] = fopen("rohit.dat", "w");
fprintf (fid, "%4.2f   %4.2f\n", table');
fclose(fid);