For \(n=3\), the Gauss nodes are \(-\sqrt {\frac {3}{5}},0,\sqrt {\frac {3}{5}}\) and the Gauss weights are \(\frac {5}{9},\frac {8}{9},\frac {5}{9}\). Write a Matlab code that verifies numerically that \(I_3 =I\) for \(f\left (x\right )=x^k\) with \(0\le k\le 5\) but \(I_3 \not = I\) for \(k=6\). Specifically, print the errors \(E_n =I_n -I\)
for each of these cases. What happens to \(E_n\) for \(k=7,8,9,10\)? How do you explain this?
nodes=[-sqrt(3/5) 0 sqrt(3/5)];
weights=[5/9 8/9 5/9];
for k=0:10
fprintf('E_%d = %0.10f\n',k,sum(nodes.^k.*weights)-integral(@(x)x.^k,-1,1))
end
When \(k\) is odd, \(f\left (x\right )=x^k\) is odd function of \(x\), therefore \(\int _{-1}^1 x^k \;\mathrm {d}x=0\), and \(E_k =\frac {5}{9}{\left (-\sqrt {\frac {3}{5}}\right )}^k +\frac {8}{9}0^k +\frac {5}{9}{\left (\sqrt {\frac {3}{5}}\right )}^k =0\), therefore \(E_k =0\);
When \(k\) is even, \(E_k \not = 0\) and we see that \(|E_k |\) increases with \(k\).
C2
Let \(A\) be the \(n\times n\) tridiagonal symmetric Jacobi matrix whose entries are all zero
apart from the numbers$$\frac{1}{\sqrt{\;1\cdot 3}},\frac{2}{\sqrt{\;3\cdot 5}},\frac{3}{\sqrt{\;5\cdot 7}},\cdots ,\frac{n-1}{\sqrt{\;\left(2n-3\right)\cdot \left(2n-1\right)}}$$in the upper-diagonal positions \(\left (1,2\right ),\left (2,3\right ),\cdots ,\left (n-1,n\right )\) of the matrix and also in the lower-diagonal
positions \(\left (2,1\right ),\left (3,2\right ),\cdots ,\left (n,n-1\right )\). Let \(\left \lbrace \lambda _j \right \rbrace \) and \(\left \lbrace v_j \right \rbrace \) be the eigenvalues and eigenvectors of \(A\) as computed e.g.
by the Matlab eig command. (In particular, each eigenvector should be
normalized so that the sum of the squares of its entries equal to 1.)
Golub and Welsh showed that the node \(x_n\) is the eigenvalue \(\lambda {\;}_j\), and the weight \(w_j\) is
twice the square of the first component of the corresponding eigenvector \(v_j\). Use
thees properties to write a Matlab function [x,w]=gaussq(n) that returns
vectors of the \(n\) Gauss quadrature nodes and weights for any \(n\). Verify that for \(n=3\),
you get the results given in (C1).
Consider the function \(f\left (x\right )=e^x \tanh \left (4\cos \left (20x\right )\right )\) for \(x\in \left \lbrack -1,1\right \rbrack \), where \(\tanh \) is as usual the hyperbolic tangent. Make a
plot of \(f\) and estimate its integral \(I\) by Gauss quadrature with \(n=1000\). Let’s call this
estimate \(I\) even though of course it is not quite exact. Now make a semilogy
plot of \(|E_n |\) against \(n\) for \(n=50,100,150,\cdots ,950,\) making sure that the plot shows clearly where the data
points lie. Comment on the shape of this curve, both its early stages and its
later stages. (As it happens, the rate of convergence is determined by the
locations of the singularities of \(f\left (x\right )\) in the complex plane.)
[x,w]=gaussq(1000);
I=sum(f(x).*w);
E=zeros(19,1);
for i=1:19
n=50*i;
[x,w]=gaussq(n);
E(i)=abs(sum(f(x).*w)-I);
end
semilogy(linspace(50,950,19),E,'-')
xlabel('n');ylabel('|E_n|');
In its early stages, the error decreases as \(n\) increases. This is because the approximation is better with more nodes. The error is smallest (\(\approx {10}^{-16}\)) when \(n=800\), because these nodes are closest to the singularities. When \(n>800\), the error increases.
C4
Determine analytically the exact integrals over [-1,1] of the functions
(You do not have to show your working.) Then make two plots with clearly
labeled curves indicating convergence of \(n\)-point Gauss quadrature for these
three integrands as a function of \(n\). One plot should be on semilogy axes, and
the other on loglog axes. Approximately how large must \(n\) be for 6-digit
accuracy in each of the three cases? What do the shapes of the curves on the
two axes indicate about convergence rates?
f=@(t)1./(1+25*t.^2);If=(2*atan(5))/5;
g=@(t)tanh(20.*(t+1/2));Ig=1;
h=@(t)abs(t);Ih=1;
Ef=zeros(19,1);Eg=Ef;Eh=Ef;
ns=linspace(3000,7000,19);i=0;
for n=ns
i=i+1;
[x,w]=gaussq(n);
Ef(i)=abs(sum(f(w))-If);
Eg(i)=abs(sum(g(w))-Ig);
Eh(i)=abs(sum(h(w))-Ih);
end
semilogy(ns,Ef,ns,Eg,ns,Eh)
xlabel('n');ylabel('|E_n|');legend('f(x)','g(x)','h(x)')
function [x,w]=gaussq(n)
array=2:n;
array=(array-1)./sqrt((2*array-3).*(2*array-1));
[V,D]=eig(diag(array,1)+diag(array,-1));
x=diag(D);
w=(2*V(1,:).^2)';
end
我修改了一下,这次上传pdf吧,比较省事
