上机作业1

2015级研究生 《计算方法》作业

2015年11月

上机作业1

数值试验3

3-1试验目的：考察不动点迭代法的局部收敛性

试验内容： 2x?ex?3?0 至少采用3种迭代法，迭代100次，考察收敛性，改变初值符号，再做迭代。分析收敛与发散的原因。

(1)迭代原理：若实数p满足p?g?p?，p称为函数g?x?的一个不动点，迭代

pn?1?g?pn?,n?0,1,...称为不动点迭代，g?x?称为迭代函数。由不动点方程建立迭代法pn?1?g?pn?,n?0,1,...，其中p0称为初值，需要预先给定。 方程2x?ex?3?0分别对应下列不同形式的不动点方程：

xxx?g(x)??x?e?3 x?g(x)?(e?3)/221??

2x?ex?3x?g3(x)?x? ?2?ex取p0??0.5,Tol?10?5,N?100，按pn?1?gi?pn?,i?1,2,3迭代，并分析收敛性。 (2)不动点迭代法代码 编制函数文件

Iteratepro.m

function y=iteratepro(x) x1=g(x); n=1;

while (norm(x1-x) >=1.0e-5) &(n<=100) x=x1;

x1=g(x);n=n+1; end x1 n

编制函数文件 不动点方程（1） function y=g(x) y=(exp(x)-3)/2;

不动点方程（2） function y=g(x) y=exp(x)-x-3; 不动点方程（3）

function y=g(x)

y=x-(2*x-exp(x)+3)/(2-exp(x)); （3）迭代结果汇总表： 近似解x* 格式（1） 格式（2） 格式（3） （4）结果分析：

（1）三种迭代格式均收敛，近似解为x*??1.3734。 （2）当设定的初值接近真实解，所用的迭代次数较小。

（3）不同迭代格式，效果不一样，迭代格式（3）实质是newton迭代法，其迭代效果明显快于其他格式。

3-2试验目的：考察Newton法求根的收敛速度

试验内容：应用Newton法求解试验3-1中的方程，并与实验3-1中收敛的迭代法进行比较，考察收敛速度。精确到0.00001。 （1）Newton迭代法原理:

在求解非线性方程f(x)?0时，它的困难在于f(x)是非线性函数，为克服这一困难，考虑它的线性展开。设当前点为xk，在xk处的Taylor展开式为

初值x=-0.5 -1.1967 -1.3734 -1.3734 迭代次数n 7 42 4 初值x=0.5 -1.3734 -1.3734 -1.3734 迭代次数n 8 41 5 f(x)?f(xk)?f'(xk)(x?xk) 令f(x)?0，解其方程得到

xk?1?xk?f(xk),(k?0,1,?) 'f(xk)

（2）程序代码 编制函数文件

newton.m

function y=newton(x0) x1=x0-fc(x0)/df(x0); n=1;

while （abs(x1-x0)>=1.0e-5）&（n<=100） x0=x1;

x1=x0-fc(x0)/df(x0);n=n+1; end x1 n

对于试验3-1编制函数文件 fc.m

function y=fc(x) y=2*x-exp(x)+3; df.m

function y=df(x) y=2-exp(x); （3）运行结果 结果汇总表： 近似解x* 格式（1） 格式（2） Newton迭代

结果分析：

（1）三者结果均收敛，解x*??1.3734。

（2）Newton迭代法迭代次数最小，收敛速度明显优于格式（1），（2）。

初值x=-0.5 -1.1967 -1.3734 -1.3734 迭代次数n 7 42 4 初值x=0.5 -1.3734 -1.3734 -1.3734 迭代次数n 8 41 5 3-3 试验目的：掌握求重根的方法

试验内容：分别用Newton法和不动点迭代法求解方程x?sinx?0，考察收敛速度，再用求重根的两种方法求方程的根，精确到0.00001。 （1）实验原理同3-1,3-2； （2）程序代码 Newton迭代法程序 函数文件 fc.m

function y=fc(x) y=x-sin(x); df.m

function y=df(x) y=1-cos(x); 程序代码

newton.m

function y=newton(x0) x1=x0-fc(x0)/df(x0); n=1;

while abs(x1-x0)>=1.0e-5 x0=x1;

x1=x0-fc(x0)/df(x0);n=n+1; end

不动点迭代法 编制函数文件

Iteratepro.m

function y=iteratepro(x) x1=g(x); n=1;

while norm(x1-x)>=1.0e-6 x=x1;

x1=g(x);n=n+1; end 函数文件 g.m

function y=g(x) y=sin(x); 运行结果

近似解 牛顿法 初值x=1 1.9449e-04 迭代次数n 21 不动点法 0.0182 k=410

(3)为了提高牛顿法求重根的收敛阶采用以下俩种方法： 方法一：

程序：function[p1,err,y]=newtonroot(f,df,p0,eps,max1) p0

for k=1:max1

p1=p0-2*（feval('f',p0)/feval('df',p0)）； err=abs(p1-p0); p0=p1;

p1,err,k,y=feval('f',p1) if(err

end 取初值p=0.5,当迭代次数k=9满足精度要求,近似解为2.4925e-05 方法二：程序：

function[p1,err,y]=newtonroot(f,df,ddf,p0,eps,max1) p0

for k=1:max1

p1=p0-(feval('f',p0)*feval('df',p0))/((feval('df',p0))^2-feval('f',p0)*feval('ddf',p0)) err=abs(p1-p0); p0=p1;

p1,err,k,y=feval('f',p1) if(err

取初值p=0.5,k =3执行结果为： p1 =1.5022e-08 err =2.2741e-08 ans = 1.5022e-08

以上四种方法中简单迭代法的收敛速度最慢 最后一种收敛速度最快。

3-4 试验目的：体验Steffensen’s method加速技巧

试验内容： 先用Newton法求解方程x?tanx?0，再用Steffensen’s method求解，比较迭代步数，精确到0.00001。

（1）简单原理：其基本思想是：对给定的初值p0(0),首先应用不动点迭代法pn+1=g(pn)计算俩步，然后用Atiken公式加速。 （2）程序代码 Newton法迭代

程序：function[p1,err,y]=newtonroot(f,df,p0,eps,max1) p0

for k=1:max1

p1=p0-feval('f',p0)/feval('df',p0); err=abs(p1-p0); p0=p1;

p1,err,k,y=feval('f',p1) if(err

定义迭代函数：存储为f.m function y=f(x) y=x-tan(x) end

定义迭代函数：存储为df.m function g=df(x) g=1-1/(cos(x))^2 end

Atiken公式加速程序：

function [n,err,p1,p2]=steffen(g,p0,tol,max1) p0

n=0,p(1)=p0; while n<=max1 for k=2:3

p(k)=feval('g',p(k-1)); end

p1=p(1)-(p(2)-p(1))^2/(p(3)-2*p(2)+p(1)); err=abs(p1-p(1)); n=n+1; p(1)=p1 n,err;

if(err

（3）结果与分析 牛顿法

取初值p0=0.5，迭代次数k=20满足精度要求. p1 =1.6018e-04 err =8.0089e-05 k =20 Ans=1.6018e-04 Atiken公式加速

取初值p0=0.5，迭代次数k=20满足精度要求。 err =8.6054e-05 n =20 p =0.1776×1.0e-03 二者迭代步数相同。从结果上看牛顿法更精确一些。

上机作业2

数值试验5-1

实验目的：熟悉Jacobi、Seidel、Sor迭代法，了解松弛因子对收敛速度的影响。

实验内容：分别用Jacobi、Seidel、Sor???0.8,1.1,1.2,1.3,1.4,1.5?迭代法求解下面的方程组，并作结果分析。

初值x(0)??0,0,0,0,0?，精度要求：x?k?1??x(k)T??10?5x(k?1)?

（1） （2）

?12.3x1?2x2?x3?3.4x4?3.7x5?4.8?13.3x1?4x2?x3?3.5x4?3.8x5?5.8?1.4x?9x?3x?2.4x?2.7x?2.3?3.4x?9x?3x?4.4x?2.3x?4.31234512345?????2.1x1?x2?8x3?2.6x4?5.8x5?2.5?4.1x1?x2?7x3?2.7x4?5.9x5?2.6?3.5x?2.1x?x?13x?4.6x?3.6?2.5x?2.4x?x?13x?5.6x?3.81234512345????2.5x1?x2?2x3?5.3x4?14.8x5?2.2??1.5x1?x2?3x3?4.3x4?14.9x5?4.2

1(1)实验原理

Jacobi迭代法求解方程组（1） 程序代码：

先建立M文件，保存文件名为jacobi.m。

function y=jacobi(a,b,x0) D=diag(diag(a)); U=-triu(a,1); L=-tril(a,-1); B=D\\(L+U); f=D\\b;

y=B*x0+f;n=1;

while norm(y-x0)>=1.0e-5 x0=y;

y=B*x0+f;n=n+1; end

y n

Seidel迭代法

D=diag(diag(A)); D0=inv(D); L=tril(A)-D; U=triu(A)-D; M=L+D; M0=inv(M); B=-M0*U;

d=M0*b;

x0=[0 0 0 0 0]'; x=B*x0+d;

while norm(x-x0)>1e-5*norm(x) x0=x; x=B*x0+d; i=i+1; end

Sor迭代式求解方程组 程序代码：

建立文件名为sor.m的M文件，具体如下： function y=sor(a,b,w,x0) D=diag(diag(a)); U=-triu(a,1); L=-tril(a,-1);

B=(D-w*L)\\((1-w)*D+w*U); f=(D-w*L)\\b*w; y=B*x0+f;n=1;

while norm(y-x0)>=1.0e-5 x0=y; y=B*x0+f; n=n+1; End y n

2程序计算结果

三种迭代法求解方程组（1）结果列表 Sor 0.8 1.1 1.2 1.3 1.4 1.5 NaN NaN NaN Jacobi — 0.3906 0.1678 0.0996 Seidel — 0.3906 0.1678 0.0996 ?x1?0.3906 0.3906 0.3906 0.3906 0.3906 ???x2?0.1678 0.1678 0.1678 0.1678 0.1678 x*??x3?0.0996 0.0996 0.0996 0.0996 0.0996 ?? ?x4??x??5?

上机作业5

7-2 实验目的 观察Lagrange插值的Runge现象， 了解若能采用合适的节点分布，则可以避免Runge现象，熟悉三次样条插值。 实验内容：对于函数f（X)=

1（-1≦x≦1）进行Lagrange插值。

1?25x^2k?，n取不同的等分数n=5,10，将区间[-1,1]等分，取等距节点。把f（x）和5,10次插值多项式的曲线画在同一张图上进行比较。再取Chebyshev节点xk=-cos

k=0,1，...，10进行Lagerange，把f（x）和Chebyshev节点的10次插值多项式的曲线画在同一张图上。 （1）原理

5.1 Lagrange插值多项式 Lagrange插值多项式的表达式：

L(x)??yili(x),i?1n?1li(x)??j?1j?in?1(x?xj)(xi?xj),i?1,2,?,n?1。

其中li(x)被称为插值基函数

（2）程序代码：

function f =Language(x,y,x0) syms t;

if(length(x) == length(y)) n = length(x); else

disp('x和y的维数不相等'); return; end

f = 0.0; for(i = 1:n) l = y(i); for(j = 1:i-1)

l = l*(t-x(j))/(x(i)-x(j));

end;

for(j = i+1:n)

l = l*(t-x(j))/(x(i)-x(j)); end;

f = f + l; simplify(f); if(i==n)

if(nargin == 3)

f = subs(f,'t',x0); else

f = collect(f);

f = vpa(f,6); end end End

（3）运行结果

（1）五等分区间

>> x=[-1 -0.6 -0.2 0.2 0.6 1];

>> y=[0.0384 0.1 0.5 0.5 0.1 0.0384]; >> f=Language(x,y)

f =1.20182*t^4 - 1.73073*t^2 + 0.567306

（2）十等分区间

>> x=[-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1];

>> y=[0.0384 0.05882 0.1 0.2 0.5 1 0.5 0.2 0.1 0.05882 0.0384]; >> f=language(x,y)

f =- 220.942*t^10 + 494.91*t^8 - 381.434*t^6 + 123.36*t^4 - 16.8552*t^2 + 1.0

（3）Chebyshev结点

>> x=[-1.0000 -0.9511 -0.8090 -0.5878 -0.3090 0 0.3090 0.5878 0.8090 0.9511 1.0000];

>> y=[0.0385 0.0424 0.0576 0.1038 0.2952 1.0000 0.2952 0.1038 0.0576 0.0424 0.0385 ]；

>> f=Language(x,y)

f =- 28.076*t^10 + 85.3475*t^8 - 96.4067*t^6 + 49.4719*t^4 - 11.2983*t^2 + 1.0

函数图像：

（1）f（x）和5,10次插值多项式曲线图像

>> t=-1:0.01:1;

>> z1=1.20182*t.^4 - 1.73073*t.^2 + 0.567306;

>> z2=- 220.942*t.^10 + 494.91*t.^8 - 381.434*t.^6 + 123.36*t.^4 - 16.8552*t.^2 + 1.0;

>> z3=(1+25*t.^2).^(-1); >> plot(t,z1,t,z2,t,z3);grid

（1） 并不是插值节点越多，插值多项式逼近函数效果就越好。 （2） 误差较大地方，是在插值区间两端点附近出现。

（2）F（x）和Chebyshev结点十次插值

>> t=-1:0.01:1;

>> z1=(1+25*t.^2).^(-1);

>> z2=- 28.076*t.^10 + 85.3475*t.^8 - 96.4067*t.^6 + 49.4719*t.^4 - 11.2983*t.^2 + 1.0; >> plot(t,z1,t,z2)

结论

利用插值多项式逼近函数，随着插值点个数的增加，其次数不断增加，由于高次多项式具有震荡特性，盲目增加节点效果反而不好。高次插值收敛没有保证，当插值区间较长时，我们宜采用分段低次插值。

上机作业6

9-1 实验目的 熟悉数值积分公式，掌握数值计算定积分的方法 实验内容：采用不同方法数值计算积分

?ln(1?x)/xdx

01（1）编写复化梯形公式和复化Simpson公式通用子程序，分别采用4,8.16,32,64等分区间计算。 （2）使用Romberg求积公式。

（3）使用高斯-勒让德求积公式（n=2,4,8）

1、试验原理：

把积分区间[a,b]分成n个子区间[xi,xi+1]，i=0，1,2,…,n，其中x0=a，xn+1=b。这样求定积分问题就分解为求和问题：

S??f(x)dx???ai?1bnxixi?1f(x)dx

当这n+1个结点为等距结点时，即xi?a?ih，其中h?(b?a)/n，i=0,1,2,…,n，复化梯形公式的形式是

hnSn??[f(xi?1)?f(xi)]

2i?1如果n还是一个偶数，则复合Simpson积分的形式是

hn/2Sn??[f(x2i?2)?4f(x2i?1)?f(x2i)]

3i?1

在区间[-1，1]上取权函数 ρ（x）=1，那么相应的正交多项式为Legendre多项式。以Legendre多项式的零点为Gauss点的求积公式为

称之为Gauss-Legendre求积公式。

2、程序代码与运行 （1）复化梯形公式

function I_n=ComTrapezoidal(a,b,n) h=(b-a)/n; for(k=0:n)

x(k+1)=a+k*h; if (x(k+1)==0)

x(k+1)=10^(-10); end end

I_1=h/2*(f(x(1))+f(x(n+1))); for (i=2:n)

F(i)=h*f(x(i)); end

I_2=sum(F); I_n=I_1+I_2;

function f=f(x) f=log(1+x)/x; 程序运行结果 n 4 8 16 32 64 1f(x)dx 0.8241 0.8231 0.8226 0.8225 0.8225 ?0

（2）复化Simpson公式程序 ComSimpson.m

function s=romberseq(fun,a,b,ep) if nargin==3 ep=1.0e-6; elseif nargin<3 error end

t1=10000; t2=-10000; n=0; m=1; h=b-a;

t(1,1)=0.5*(b-a)*(feval(fun,a)+feval(fun,b)); while abs(t2-t1)>=ep

area=0.0; n=n+1; h=h/2; for i=1:m

area=area+feval(fun,h*(2*i-1)+a); end

t(n+1,1)=0.5*t(n,1)+area*h; m=2*m; if n>4

for j=1:3

for i=1:n-j

t(i,j+1)=(4^(j)*t(i+1,j)-t(i,j))/(4^(j)-1); end end

t1=t(n-4,4); t2=t(n-3,4); end end disp; s=t2;

建立脚本文件

exf=inline('(log(1+x))./x'); fplot(exf,[0,1]) comsimp(exf,0,1,4)

运行结果: n 4 8 16 32 64 1f(x)dx 0.8225 0.8225 0.8225 0.8225 0.8225 ?0

（3）Romberg求积公式程序

a=0;b=1;eps=0.00001;R=zeros(4,4); h=b-a;err=1;J=0;

R(1,1)=h*(tx(a)+tx(b))/2; while (err>eps) J=J+1; h=h/2;

x=((a+h):2*h:(b-h)); [m,n]=size(x); A=zeros(n,1); for i=1:n

A(i,1)=tx(x(m,i)); end

R(J+1,1)=R(J,1)/2+h*sum(A); for K=1:min(3,J)

R(J+1,K+1)=R(J+1,K)+(R(J+1,K)-R(J,K))/(4^K-1); end

if (J>3)

err=abs(R(J+1,4)-R(J,4)); end end

quad=R(J,4); 步长 1 T 0.846573590279973 0.5 0.828751903248151 0.25 0.824058098916759 0.125 0.822866129037621 0.0625 0.822566892025848 0.822811340904210 0.822493497472962 0.822468805744575 0.822467146355257 0.822472307910879 0.822467159629349 0.822467035729302 0.822467077910594 0.822467033762635 S C R 积分近似值为0.822467033762635。

（4）高斯-勒让德求积公式程序 guass1.m

function s=guass1(a,b,n) h=(b-a)/n; s=0.0;

for m=0:(1*n/2-1)

s=s+h*(guassf(a+h*((1-1/sqrt(3))+2*m))+guassf(a+h*((1+1/sqrt(3))+2*m)));

end

guassf.m

function y=guassf(x) y =log(1+x)/x;

高斯-勒让德求积公式结果汇总： 2 n 0.8222 1?f(x)dx 04 0.8224 8 0.8225 结论

分析并比较得到的数据可以看出，当k越来越大时，数值积分的结果越来越靠近原函数积分实际结果，并且复合Simpson积分的结果更迅速地靠近原函数积分实际结果，这是有原因的，从两种方法的误差项即可看出。

复合梯形积分的误差项是?是?1(b?a)h4f180(4)1(b?a)h2f''(?)，Simpson积分的误差项12(?)，??(a,b)，当h趋于零时，显然Simpson积分的

误差项更快地趋于零，实验结果复符合这一结论。

本文来源：https://www.bwwdw.com/article/2rlg.html