科学工程计算Word格式.docx

科学工程计算Word格式.docx

《科学工程计算Word格式.docx》由会员分享，可在线阅读，更多相关《科学工程计算Word格式.docx（22页珍藏版）》请在冰豆网上搜索。

【题目】使用牛顿迭代法和斯蒂芬森（Steffensen）加速法求解x^3-5x^2+6x-1=0全部根，（注：

根全部在区间【0,4】中）要求精确到10^（-6） 

并给出以[-5,5]区间的100等分点（即-5，-4.9,-4.8,.....,4.8,4.9,5）为初始点，收敛到不同解的示意图

算法说明：

牛顿法递推公式为：

（1.1）

取-5，-4.9,-4.8,.....,4.8,4.9,5为初始点，分别代入式1.1，求出满足条件的解（由于题目中给出其根全部在区间【0,4】中，此方法求出的解即为x^3-5x^2+6x-1=0全部根）

源代码：

tol=10^（-6）;

%误差限

N=2000;

%最大迭代步骤

i=1;

forp0=-5:

0.1:

5;

fork=1:

N

p1=p0-（p0^3-5*p0^2+6*p0-1）./（3*p0^2-10*p0+6）;

ifabs（p1-p0）<

tol

break

end

p0=p1;

end

a（i）=p1;

i=i+1;

end

p0=-5:

5

c0=find（abs（a-0.1981）<

0.1）,c1=find（abs（a-1.555）<

0.1）,c2=find（abs（a-3.2470）<

0.001）

subplot（111）

plot（p0（c0）,0,'

b-*'

p0（c1）,0,'

r-d'

p0（c2）,0,'

c-o'

）;

grid

axis（[04-1010]）;

输出结果：

1.2斯蒂芬森（Steffensen）加速法

Steffensen加速法递推公式为：

，

且有

n=0;

p

（1）=p0;

whilen<

=N

fork=1:

2

p（k+1）=p（k）-（p（k）^3-5*p（k）^2+6*p（k）-1）/（3*p（k）^2-10*p（k）+6）;

%牛顿迭代法

p1=p

（1）-（p

（2）-p

（1））^2/（p（3）-2*p

（2）+p

（1））;

%Aitken加速

ifabs（p1^3-5*p1^2+6*p1-1）<

tol

break

n=n+1;

p

（1）=p1;

线性方程组数值解法程序设计与验证

【实验目的及意义】

1）掌握Gauss消去法及Gauss列主元消去法，能用这两种方法求解方程组。

2）掌握矩阵和向量范数的定义，了解它们的性质。

3）掌握Jacobi迭代法，能应用Jacobi迭代法求解方程组；

4）掌握Seidel迭代法，能应用Seidel迭代法求解方程组。

5）能用矩阵的范数判断迭代法的收敛性。

【实验要求】

设计和验证Gauss按比例列选主元消去法和Jacobi迭代法、Seidel迭代法算法的程序

误差不超过

.

2.1Gauss按比例列选主元消去法

为了克服Gauss消元法的两点不足，在每步消元时先选择绝对值较大的数作为乘数mik的分母，控制舍入误差的扩大。

列主元消元法是一种局部主元消元法，按比例列选主元消去法是在列主元消元法的基础上给出的一种改进，主要作用是消除方程中比例因子对选主元的影响。

A=[17.031-0.615-2.9111.007-1.0060.000;

-1.00034.211-1.000-2.1000.300-1.700;

0.0000.50013.000-0.5001.000-1.500;

4.5013.110-3.907-61.70512.1708.999;

0.101-8.012-0.017-0.9104.6180.100;

1.000234.5521.803];

b=[0.230;

-52.322;

54;

240.236;

29.304;

-117.818];

B=[A,b]

Tol=10^（-6）;

[n,m]=size（B）;

ifn~=m-1

'

error'

return

fori=1:

n

W（i）=max（B（i,:

））;

n-1

K=B（k:

n,k）'

./W（k:

n）;

[l,m]=max（K）;

ifm~=1

T=B（k,:

B（k,:

）=B（k+m-1,:

B（k+m-1,:

）=T;

S=W（k）;

W（k）=W（k+m-1）;

W（k+m-1）=S;

fori=k+1:

1:

B（i,:

）=B（i,:

）-B（k,:

）*B（i,k）/B（k,k）;

ifB（n,n）==0

elseX（n）=B（n,n+1）/B（n,n）;

fork=n-1:

-1:

1

X（k）=（B（k,n+1）-dot（B（k,k+1:

n）,X（k+1:

n）））/B（k,k）;

ifabs（X（k）-X（k+1））<

Tol

end

disp（X）;

2.2Jacobi迭代法

上两式分别为Jacobi不动点方程、Jacobi迭代格式。

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;

whilenorm（y-x0）>

=1.0e-6

x0=y;

y=B*x0+f;

n=n+1;

y

2.3Seidel迭代法

x0=[0;

0;

0];

N=100;

[n,n]=size（A）;

x=x0;

k=1;

whilek<

fori=1:

x（i）=（b（i）-A（i,:

）*x+x（i）*A（i,i））/A（i,i）;

ifabs（x-x0）<

break;

k=k+1;

x0=x;

disp（x）

曲线拟合

x=0,1,2,3,4,5,6,7,8,910

y=7.579,8.454,6.819,4.620,3.161,3.108,4.491,6.698,8.482,7.956,2.594

使用以下方法求0到10的100等分点（即0.1,0.2,0.3,........9.8,9.9）的插值或拟合值

（1）牛顿插值（11点全用）

（2）任何点选择离其最近的四个节点做拉格朗日差值（如点3.4使用节点2,3,4,5进行插值，8.8使用节点7,8,9,10插值

（3）三次拟合曲线

（4）四次拟合曲线

（5）上述四种插值或拟合值在同一图上画图比较

3.1牛顿插值

x=[012345678910];

y=[7.5798.4546.8194.6203.1613.1084.4916.6988.4827.9562.594];

n=10;

forj=1:

fori=n+1:

j+1

y（i）=（y（i）-y（i-1））/（x（i）-x（i-j））;

y1=y（n+1）;

x1=0.1;

99

forj=n:

1

y1=y（j）+（x1-x（j））*y1;

x2（k）=x1;

y2（k）=y1;

x1=x1+0.1;

3.2拉格朗日插值

首先将x=1:

10的100的等分，得101个值，每一个值取整，求其最近的4个数来做插值，这样每一个值都得到4个数为一组的插值，然后利于拉格朗日插值公式：

Ln（x）=

将每组数的4个值代入得到各自的插值基函数，并将其对应的x1=0:

10中的值代入到所求基函数中得到对应的插值，及存入到数组y1中。

x1=0:

10;

101

y1（k）=0;

i=floor（x1（k））;

ifi<

form=0:

3

t=1;

forj=0:

ifj~=m

t=t*（x1（k）-x（j+1））/（x（m+1）-x（j+1））;

y1（k）=y1（k）+t*y（m+1）;

elseifi<

9

form=i-1:

i+2

forj=i-1:

else

form=7:

10

forj=7:

3.3三次拟合曲线

取

，然后根据公式

分别求的法矩阵A4x4,由于根据

求的d矩阵，这样解得法方程A*a=d,解得a矩阵，及求的

所以三次拟合曲线

将x1=0:

10的101个点代入到三次拟合曲线多项式中解得101个拟合值，并存入到数组y1中。

A=zeros（4,4）;

d=zeros（4,1）;

A（1,1）=11;

A（1,2）=sum（x）;

A（2,1）=sum（x）;

A（1,3）=sum（x.^2）;

A（3,1）=sum（x.^2）;

A（2,2）=sum（x.^2）;

A（1,4）=sum（x.^3）;

A（4,1）=sum（x.^3）;

A（2,3）=sum（x.^3）;

A（3,2）=sum（x.^3）;

A（2,4）=sum（x.^4）;

A（4,2）=sum（x.^4）;

A（3,3）=sum（x.^4）

A（3,4）=sum（x.^5）;

A（4,3）=sum（x.^5）;

A（4,4）=sum（x.^6）;

d

（1）=sum（y）;

d

（2）=sum（x.*y）;

d（3）=sum（x.^2.*y）;

d（4）=sum（x.^3.*y）;

a=A\d;

y1（k）=a（4）;

fori=3:

y1（k）=a（i）+x1（k）*y1（k）;

3.4四次拟合曲线 

分别求的法矩阵A5x5,由于根据

所以四次拟合曲线

10的101个点代入到四次拟合曲线多项式中解得101个拟合值，并存入到数组y1中。

A=zeros（5,5）;

d=zeros（5,1）;

A（1,4）=sum（x.^3）;

A（2,3）=sum（x.^3）;

A（2,4）=sum（x.^4）;

A（4,2）=sum（x.^4）;

A（1,5）=sum（x.^4）;

A（5,1）=sum（x.^4）;

A（3,3）=sum（x.^4）;

A（4,3）=sum（x.^5）;

A（2,5）=sum（x.^5）;

A（5,2）=sum（x.^5）;

A（4,4）=sum（x.^6）;

A（5,3）=sum（x.^6）;

A（3,5）=sum（x.^6）;

A（5,4）=sum（x.^7）;

A（4,5）=sum（x.^7）;

A（5,5）=sum（x.^8）;

d

（2）=sum（x.*y）;

d（3）=sum（x.^2.*y）;

d（4）=sum（x.^3.*y）;

d（5）=sum（x.^4.*y）;

y1（k）=a（5）;

fori=4:

3.5四种插值拟合值图形比较

plot（x1,newton2（）,'

b:

*'

x1,langange2（）,'

g-.+'

x1,lineUntitled3（）,'

r-x'

x1,lineUntitled4（）,'

k--d'

legend（'

牛顿插值'

'

拉格朗日插值'

三次曲线拟合'

四次曲线拟合'

Location'

North'

几种求积公式积分

【题目】

（1）使用第三次作业的四次拟合函数的拟合值

分别使用复化梯形公式和复化simpson公式求（一共101点的数值）的积分近似值

（2）编程实现书上第九章第6题

p=polyfit（x,y,4）;

a=0;

b=10;

fa=polyval（p,a）;

fb=polyval（p,b）;

f0=fa+fb;

h=（b-a）/100;

y1=a+h:

h:

b-h;

f1=sum（polyval（p,y1））;

T=h/2*（f0+2*f1）

2*h:

y2=a+2*h:

b-2*h;

f2=sum（polyval（p,y2））;

T=h/3*（f0+4*f1+2*f2）

4.3龙贝格积分法计算

symsx;

b=4;

h=b-a;

eps==10^（-5）;

err=1;

J=0;

R=zeros（4,4）;

f=exp（x）;

fa=subs（f,'

x'

a）;

fb=subs（f,'

b）;

R（1,1）=h*（fa+fb）/2;

while（err>

eps）

J=J+1;

h=h/2;

x=a+h:

R（J+1,1）=R（J,1）/2+h*sum（subs（f,'

x））;

min（3,J）

R（J+1,k+1）=R（J+1,k）+（R（J+1,k）-R（J,k））/（4^k-1）;

if（J>

3）

err=abs（R（J+1,4）-R（J,4））;

quad=R（J,4）

常微分方程初值问题的数值解法

分别使用改进euler方法和四级四阶RK方法求解10－3题

取步长h＝0.01，计算到y（0.5）

5.1改进euler方法

h=0.01;

n=50;

x=0;

y=0;

f1=-x^2+x-y;

y1=y+h*f1;

f2=-（x+h）^2+x+h-y1;

y2=y+h*f2;

y=（y1+y2）/2;

a（i）=y;

x=x+h;

disp（a）;

5.2四级四阶RK方法

对于微分方程，经典龙格库塔法的计算公式为：

k1=h*f1;

f2=-（x+h/2）^2+x+h/2-y-k1/2;

k2=h*f2;

f3=-（x+h/2）^2+x+h/2-y-k2/2;

k3=h*f3;

f4=-（x+h）^2+x+h-y-k3;

k4=h*f4;

y=y+（k1+2*k2+2*k3+k4）/6;

b（i）=y;

disp（b）;