Section 2 Problem 16 in Greenspan and Casulli
The answer given in Greenspan and Casulli
for problem 16 on page 72 is different than the answer
I obtain using a Maple worksheet to solve the normal equations.
The last graph compares my answer (in green) and
the book's answer (in blue) with the data.
|\^/| Maple V Release 5 (University of Nevada, Reno)
._|\| |/|_. Copyright (c) 1981-1997 by Waterloo Maple Inc. All rights
\ MAPLE / reserved. Maple and Maple V are registered trademarks of
<____ ____> Waterloo Maple Inc.
| Type ? for help.
> with(linalg);
Warning, new definition for norm
Warning, new definition for trace
[BlockDiagonal, GramSchmidt, JordanBlock, LUdecomp, QRdecomp, Wronskian, addcol,
addrow, adj, adjoint, angle, augment, backsub, band, basis, bezout,
blockmatrix, charmat, charpoly, cholesky, col, coldim, colspace, colspan,
companion, concat, cond, copyinto, crossprod, curl, definite, delcols,
delrows, det, diag, diverge, dotprod, eigenvals, eigenvalues, eigenvectors,
eigenvects, entermatrix, equal, exponential, extend, ffgausselim, fibonacci,
forwardsub, frobenius, gausselim, gaussjord, geneqns, genmatrix, grad,
hadamard, hermite, hessian, hilbert, htranspose, ihermite, indexfunc,
innerprod, intbasis, inverse, ismith, issimilar, iszero, jacobian, jordan,
kernel, laplacian, leastsqrs, linsolve, matadd, matrix, minor, minpoly,
mulcol, mulrow, multiply, norm, normalize, nullspace, orthog, permanent,
pivot, potential, randmatrix, randvector, rank, ratform, row, rowdim,
rowspace, rowspan, rref, scalarmul, singularvals, smith, stackmatrix,
submatrix, subvector, sumbasis, swapcol, swaprow, sylvester, toeplitz,
trace, transpose, vandermonde, vecpotent, vectdim, vector, wronskian]
> Digits:=10;
Digits := 10
> alpha:=[.02802,.03040,.03247,.02850,.02900,.02497,
> .02447,.02856,.02566];
alpha :=
[.02802, .03040, .03247, .02850, .02900, .02497, .02447, .02856, .02566]
> R:=[2910.0,1860.0,1380.0,1100.0,910.0,780.0,680.0,600.0,540.0];
R := [2910.0, 1860.0, 1380.0, 1100.0, 910.0, 780.0, 680.0, 600.0, 540.0]
> A:=matrix([seq([1,R[i],R[i]^2]*alpha[i],i=1..9)]);
[.02802 81.538200 237276.1620]
[.03040 56.544000 105171.8400]
[.03247 44.808600 61835.86800]
[.02850 31.350000 34485.00000]
A := [.02900 26.390000 24014.90000]
[.02497 19.476600 15191.74800]
[.02447 16.639600 11314.92800]
[.02856 17.136000 10281.60000]
[.02566 13.856400 7482.456000]
> At:=transpose(A);
At :=
[.02802 , .03040 , .03247 , .02850 , .02900 , .02497 , .02447 , .02856 ,
.02566]
[81.538200 , 56.544000 , 44.808600 , 31.350000 , 26.390000 , 19.476600 ,
16.639600 , 17.136000 , 13.856400]
[237276.1620 , 105171.8400 , 61835.86800 , 34485.00000 , 24014.90000 ,
15191.74800 , 11314.92800 , 10281.60000 , 7482.456000]
> N:=multiply(At,A);
[.0071132223 8.855819304 14674.62379 ]
[ ]
[ 8 ]
N := [8.855819304 14674.62379 .3054356919 10 ]
[ ]
[ 8 11]
[14674.62379 .3054356919 10 .7347121346 10 ]
> Y:=[seq(log(alpha[i])*alpha[i],i=1..9)];
Y := [-.1001669254, -.1061967052, -.1112889344, -.1013987590, -.1026733240,
-.09214130197, -.09079122213, -.1015521669, -.09398801047]
> b:=multiply(At,Y);
b := [-.02535049158, -31.39518107, -51954.66367]
> Ni:=inverse(N);
[ 2921.158092 -4.071088893 .001108987934 ]
[ ]
[ -5]
Ni := [-4.071088893 .006179512782 -.1755826157 10 ]
[ ]
[ -5 -9 ]
[.001108987934 -.1755826157 10 .5220441467 10 ]
> c:=multiply(Ni,b);
[ -6]
c := [-3.85731585, .00042053945, -.11153721 10 ]
> T1:=solve({c[1]=log(A0*exp(-B0*C0^2)),c[2]=2.0*B0*C0,c[3]=-B0},
> {A0,B0,C0});
-6
T1 := {B0 = .1115372100 10 , C0 = 1885.197998, A0 = .03140099228}
> P:=[seq([R[i],alpha[i]],i=1..9)];
P := [[2910.0, .02802], [1860.0, .03040], [1380.0, .03247], [1100.0, .02850],
[910.0, .02900], [780.0, .02497], [680.0, .02447], [600.0, .02856],
[540.0, .02566]]
> f:=unapply(subs(T1,A0*exp(-B0*(x-C0)^2)),x);
-6 2
f := x -> .03140099228 exp(-.1115372100 10 (x - 1885.197998) )
> g:=x->3.05*10^(-2)*exp(-2.63*10^(-8)*(x-2880)^2);
-7 2
g := x -> .03050000000 exp(-.2630000000 10 (x - 2880) )
> plot([P,f(x),g(x)],x=500..3000);
> Ef:=sum((alpha[i]-f(R[i]))^2,i=1..9);
Ef := .00002291952637
> Eg:=sum((alpha[i]-g(R[i]))^2,i=1..9);
Eg := .00003772393277
> quit
bytes used=2063172, alloc=1441528, time=3.38