《现代数值计算》Matlab程序整理

7
2.将积分区间 n 等分，用复合梯形求积公式计算定积分 （画出平面上 log(n)-log(Error)图） ． 解： >> edit ksf7.m function v=ksf7(x) v=sqrt(1+x.^2); >> I=quad('ksf7',1,3); >> n=1:100; >> for i=1:100 x(i)=tixingjifen('ksf7',1,3,i); error(i)=abs(I-x(i)); end >> plot(log10(n),log10(error(n)))
4
第 3 章 多项式插值与样条插值
1.拉格朗日插值
function yh=lagrange(x,y,xh) n=length(y); yh=zeros(size(xh)); for k=1:n, pt=ones(size(xh)); for j=[1:k-1 k+1:n], pt=pt.*(xh-x(j))/(x(k)-x(j)); yh=yh+y(k)*pt; end end >> x=[11,12]; >> y=[2.3979,2.4829]; >> xh=[11.75]; >> yh=lagrange(x,y,xh) yh = 2.4616
2
2. LU 分解 function [L,U]=myLU(A) %实现对矩阵 A 的 LU 分解，L 为下三角矩阵 [n,n]=size(A); L=zeros(n,n); U=zeros(n,n); for i=1:n L(i,i)=1; end for k=1:n for j=k:n U(k,j)=A(k,j)-sum(L(k,1:k-1)*U(1:k-1,j)); end for i=k+1:n L(i,k)=(A(i,k)-sum(L(i,1:k-1)*U(1:k-1,k)))/U(k,k); end end 3.追赶法求解线性方程组，n=101.
1
1.2 写一个求解线性代数方程 A=[1,13,-2,-34;2,6,-7,-10;-10,-1,5,9;-3,-5,0,15];b=[13,-22,14,-36]'; 的列主元素高斯消去法程序 解： function x=gaussel(A,b) n=length(b); for i=1:n-1 [r,p]=max(abs(A(i:n,i))); %r 为最大值，p 为所在行 p=p+i-1; if p>i t=A(i,:);A(i,:)=A(p,:);A(p,:)=t; t=b(i);b(i)=b(p);b(p)=t; end for j=i+1:n r=-A(j,i)/A(i,i); A(j,:)=A(j,:)+r*A(i,:); b(j)= b(j)+r*b(i); end end x=zeros(n,1); x(n)=b(n)/A(n,n); for k=(n-1):-1:1 x(k)=(b(k)-A(k,k+1:n)*x(k+1:n))/A(k,k); end >> A=[1,13,-2,-34;2,6,-7,-10;-10,-1,5,9;-3,-5,0,15]; >> b=[13,-22,14,-36]'; >> x=gaussel(A,b) x= 14.0669 29.9869 17.9950 10.4090
6
第 5 章 数值积分与数值微分
1.复合梯形公式计算下面积分的近似值，区间等分 n 20 。 复合辛普森公式、复合中点公式计算下面积分的近似值，区间等分 n 10 。
1 0
1 x
2
1
dx
解： >> edit tixingjifen.m function s=tixingjifen(f,a,b,n) x=linspace(a,b,(n+1)); y=zeros(1,length(x)); y=feval(f,x); h=(b-a)./n; s=0.5*h*(y(1)+2*sum(y(2:n))+y(n+1)); end >> edit simpson.m function I=simpson(f,a,b,n) h=(b-a)/n; x=linspace(a,b,2*n+1); y=feval(f,x); I=(h/6)*(y(1)+2*sum(y(3:2:2*n-1))+4*sum(y(2:2:2*n))+y(2*n+1)); >> edit fmid.m function I= fmid(f,a,b,n) h=(b-a)./n; x=linspace(a+h/2,b-h/2,n); y=feval(f,x); I=h*sum(y); end >> edit ksf4.m function v=ksf4(x) v=1./(x.^2+1); >> tixingjifen('ksf4',0,1,20) ans = 0.7853 >>simpson('ksf4',0,1,10) ans = 0.7854
3
4.乔列斯基 Choleskey 分解 function x=choleskey(A,b) [N,N]=size(A); x=zeros(N,1); y=zeros(N,1); for i=1:N A(i,i)=sqrt(A(i,i)-A(i,1:i-1)*A(i,1:i-1)'); for j=i+1:N A(j,i)=(A(j,i)-A(j,1:i-1)*A(i,1:i-1)')/A(i,i); end end A for j=1:N y(j)=(b(j)-A(j,1:j-1)*y(1:j-1))/A(j,j); end y A=A'; for k=N:-1:1 x(k)=(y(k)-A(k,k+1:N)*x(k+1:N))/A(k,k); end
第 2 章 线性方程组的直接解法 1.1 写一个求解线性代数方程 x b 的列主元素高斯消去法程序，在程序入口输入 n 的
值，其中 ij
1 ， b j ln j ， 1 i, j n i j 1
解： function x=gauss(n) for i=1:n for j=1:n A(i,j)=1/(i+j-1); end b(i,1)=log(i); end %Gauss Elimination for i=1:n-1 [r,p]=max(abs(A(1:n,i))); p=p+i-1; A([p,i],:)=A([i,p],:); for j=i+1:n r=-A(j,i)/A(i,i); A(j,:)=A(j,:)+r*A(i,:); b(j)= b(j)+r*b(i); end end x=zeros(n,1); for k=n:-1:1 x(k)=(b(k)-A(k,k+1:n)*x(k+1:n))/A(k,k); end
y(1 cx2 ) ax b ，进一步变形为 y ax b cx 2 y ，则基函数是 x ，1， x 2 y 。
>> x=[-0.931 -0.586 -0.362 -0.213 0.008 0.544 0.628 0.995]'; >> y=[0.356 0.606 0.687 0.802 0.823 0.801 0.718 0.625]'; >> A=[x ones(8,1) -x.^2.*y]; >> z=A\y; >> a=z(1); b=z(2); c=z(3); >> xh=[-1:0.1:1]'; >> yh=(a.*xh+b)./(1+c.*xh.^2); >> plot(x,y,'r+',xh,yh,'b*') 2.最小二乘法拟合形如 f ( x ) 拟合效果。 次序 年份 人口（亿） 解： 1 1953 5.82 2 1964 6.95 3 1982 10.08 4 1990 11.34 5 2000 12.66
2 3.函数 polyfit 最小二乘多项式拟合 p ( x) a0 a1 x a2 x
>> x=[0,0.25,0.50,0.75,1.00]; >> y=[1.00,1.284,1.6487,2.1170,2.7183]; >> p=polyfit(x,y,2)；%幂次从高到低排列，2 为多项式的阶数 >> a0=p(3); a1=p(2); a2=p(1);
用 Jacobi 迭代方法对下面方程组求解，取初始向量 x(0) (3, 2, 1)T 。
2 4 4 x1 2 3 3 3 x 3 2 x 2 3 4 4 2
解： >>edit Jacobi.m function[x iter]=Jacobi(A,x0,b,tol) D=diag(diag(A)); L=D-tril(A); U=D-triu(A); x=x0; iter=0; while(norm(A*x-b)/norm(b)>tol) x=D\((L+U)*x0+b); x0=x; iter=iter+1; end >> A=[2 4 -4;3 3 3;4 4 2]; >> b=[2 -3 -2]'; >> x0=[3 2 -1]'; >> [x,iter]=Jacobi(A,x0,b,1e-4) x= 1 -1 -1 iter = 3

3 1
1 x 2 dx ,比较不同 n 值时的误差
8
第 6 章 线性方程组的迭代解法
1. 雅可比迭代法 迭代格式： xi
( k 1) n 1 (bi aij x(jk ) ) ；迭代矩阵 BJ =D 1 ( L U ) I D 1 A aii j 1, j i
Leabharlann Baidu

5
第 4 章 函数逼近 ax b 1.下面的数据表近似地满足函数 y ，请适当变换成为线性最小二乘问题，编程求 1 cx 2
最好的系数 a, b, c ，并在同一个图上画出所有数据和函数图像.
xi yi
解：
0.931 0.586 0.362 0.213 0.008 0.544 0.628 0.995 0.356 0.606 0.687 0.802 0.823 0.801 0.718 0.625
解： function x=tridiagsolver(a,b,c,d) n=length(b); l(1)=b(1); y(1)=d(1)/l(1); u(1)=c(1)/l(1); for i=2:n-1, l(i)=b(i)-a(i-1)*u(i-1); y(i)=(d(i)-y(i-1)*a(i-1))/l(i); u(i)=c(i)/l(i); end l(n)=b(n)-a(n-1)*u(n-1); y(n)=(d(n)-y(n-1)*a(n-1))/l(n); x(n)=y(n); for i=n-1:-1:1 x(i)=y(i)-u(i)*x(i+1); end >> n=101; >> x=tridiagsolver(ones(1,n-1),12*ones(1,n),ones(1,n-1),[11,10*ones(1,n-2),11])
12 1 0 ... 0 x1 11 1 12 1 ... 0 x2 10 0 1 12 ... 0 x3 10 ... ... ... ... 1 x4 ... 0 0 ... 1 12 x 11 5
a bx 的函数，拟合给出的中国人口数据，并通过图形展示 1 cx
f ( x)(1 cx) a bx ， xf ( x) 。 进一步变形为 f ( x) a bx cxf ( x) ， 则基函数是 1，x ，
>> x=[1953,1964,1982,1990,2000]'; >> y=[5.82,6.95,10.08,11.34,12.66]'; >> A=[ones(5,1) x -x.*y]; >> z=A\y; >> a=z(1); b=z(2); c=z(3)； >> v=linspace(1953,2000,100); >> plot(x,y,'b-+',v,(a+b*v)/(1+c*v),'k-'); %blue,black