数值分析上机题参考答案.docx

如有帮助欢迎下载支持

数值分析上机题

姓名：陈作添

学号： 040816

习题 1

20．（上机题）舍入误差与有效数

N

1

1 3 1 1

设

S N

，其精确值为 。

2

2 2 N

N 1

j 2

j

1

（1）编制按从大到小的顺序

1

1

1 ，计算 S 的通用程序。

S N

1 32

1

N 2

1 N

2

2

（2）编制按从小到大的顺序

1

1

1

，计算 S 的通用程序。

S N

1

(N 1)2 1

22 1

N

N 2

（3）按两种顺序分别计算

S 102 ， S 104 ， S 106 ，并指出有效位数。 （编制程序时用单精度）

（4）通过本上机题，你明白了什么？

按从大到小的顺序计算 S N 的通用程序为： 按从小到大的顺序计算 S N 的通用程序为：

#include #include float sum(float N) float sum(float N) {

{

float j,s,sum=0; float j,s,sum=0; for(j=2;j<=N;j++) for(j=N;j>=2;j--) {

{

s=1/(j*j-1); s=1/(j*j-1); sum+=s;

sum+=s;

}

}

return sum;

return sum;

}

}

从大到小的顺序的值

从小到大的顺序的值

精确值

有效位数

从大到小

从小到大

0.740049

0.74005

0.740049

6

5 S 102

0.749852

0.7499

0.7499

4 4

S 104

0.749852

0.749999

0.749999

3

6

S 106

通过本上机题， 看出按两种不同的顺序计算的结果是不相同的，

按从大到小的顺序计算

的值与精确值有较大的误差， 而按从小到大的顺序计算的值与精确值吻合。 从大到小的顺序

计算得到的结果的有效位数少。

计算机在进行数值计算时会出现“大数吃小数”的现象，导

致计算结果的精度有所降低，

我们在计算机中进行同号数的加法时，

采用绝对值较小者先加

的算法，其结果的相对误差较小。

习题 2

20．（上机题） Newton 迭代法

（1）给定初值x 0及容许误差，编制 Newton 法解方程f ( x)0 根的通用程序。

（2）给定方程f ( x) x3/ 3x0 ，易知其有三个根x 3 ，x2 0， x3 。

13 1．由 Newton 方法的局部收敛性可知存在0 ，当x0(,) 时，Newton迭代序列收敛于根 x2。试确定尽可能大的。

2．试取若干初始值，观察当x0 (, 1)，(1, )，(,) ， (,1)，(1, ) 时Newton 序列是否收敛以及收敛于哪一个根。

（3）通过本上机题，你明白了什么？

解：（ 1）编制的通用程序：

#include{

#includefloat x0,x1,a;

#define eps0.000001/给定容许误差int k=0;

float f(float x)//定义函数 f(x)cout<<" 请输入初值 x0:";

{cin>>x0;

float f;do

f=x*x*x/3-x;//f(x) 的表达式 ;{

return(f);a=-f(x0)/df(x0);

}x1=x0+a;

float df(float x)//定义函数 df(x) ，k++;

计算 f(x) 的导函数x0=x1;

{}

float df;while(fabs(a)>eps);

df=x*x-1;//f(x) 导函数的表达式 ;cout<

return (df);//输出迭代的次数和根值

}}

void main(void)

（ 2）计算迭代序列收敛于根x2的尽可能大的的函数为：

#includereturn(f);

#include}

void delay(int n)// 定义延时函数float df(float x)// 定义函数 df(x) ，计算 f(x) {for(n=10000;n>0;n--);}的导函数

#define eps 0.000001{

float f(float x)//定义函数 f(x)float df;

{df=x*x-1;//f(x) 导函数的表达式 ;

float f;return (df);

f=x*x*x/3-x;//f(x) 的表达式 ;}

int judgement(float z)return 0;

{}

int count=5;void main(void)

float x0,x1,type,type1;{

x0=z;float delta=0;

while(count-->0)int flag=1;

{while(flag==1)

x1=x0-f(x0)/df(x0);{

type=fabs(x1);cout<<" 方程的根为 :"<<'\n';

type1=fabs(x1-x0); // 调试值用delta+=eps;

cout<<"count="<

"<

if(fabs(x1-x0)

x0=x1;cout<

delay(30000);//调试值用}

}

当步长为 0.001 时，程序计算出的δ的为δ =0.774 ，即在区间（ -0.774， 0.774）内迭代序列收敛于 0。

对于不同得初始值收敛于不同的根, x0在（ -∞， -1）内收敛于x1*，在（ -0.774， 0.774）内

收敛于x2，在（ 1，+∞）内收敛于x3，但在内（ 0.774，1）和（－ 1，0.774）均可能收敛于

x1*和 x3。 x1*， x2， x3分别为方程的精确解。

分析：对于不同的初值，迭代序列会收敛于不同的根，所以在某个区间内求根对于初值的选取

有很大的关系。产生上述结果的原因是区间不满足大范围收敛的条件。

习题 3

35．（上机题）列主元三角分解法对于某电路的分析，归结为求解线性方程组RI=V。

（1）编制解 n 阶线性方程组 Ax=b 的列主元三角分解法的通用程序；

（2）用所编制的程序解线性方程组RI=V，并打印出解向量，保留五位有效数；

（3）本编程之中，你提高了哪些编程能力？

程序为：

#includecin>>n;

#includecout<<" 输入数组 a:"<

void main(void)for(i=1;i<=n;i++)

{for(j=1;j<=(n+1);j++)

int i,j,n,k,q;cin>>a[i][j];//给矩阵 a 赋值

float a[10][11],s[10],s1[10];for(i=1;i<=n;i++)

cout<<" 请输入 n 的值 :";{

for(j=1;j<=(n+1);j++)

cout<

cout<<'\n';

}//输出数组a

cout<<"'''''''''''''''''''''''''"<<'\n';

//进行第一行和第一列元素的求取

'''''''''''''''''''''''''//

int t=1;

for(i=1;i<=n;i++)

{s[i]=a[i][1];}

float max=fabs(s[1]);

for(i=2;i<=n;i++)

if(fabs(s[i])>max)

{

max=fabs(s[i]);

t=i;

}

for(j=1;j<=(n+1);j++)

{

float b=a[1][j];

a[1][j]=a[t][j];

a[t][j]=b;

}//进行第一列主元互换for(i=2;i<=n;i++)

a[i][1]=a[i][1]/max; //第一列除以 a[1][1] for(i=1;i<=n;i++)

{

for(j=1;j<=(n+1);j++)

cout<

cout<<'\n';

}

//输出进行第一步变换的数组 a

cout<<"'''''''''''''''''''''''''"<<'\n';

//进行第 k 步分解 '''''''''''''''''''''''''''''''''''''''''// for(k=2;k<=n;k++)

{

for(i=k;i<=n;i++)

{

float sum=0;

for(q=1;q

sum+=a[i][q]*a[q][k];

s1[i]=a[i][k]-sum;

}

int l=k;float m=fabs(s1[k]);

for(i=k;i<=n;i++)

// 比较第 k 步分解的第k 列值的大小{

if(fabs(s1[i])>m)

{

m=fabs(s1[i]);

l=i;// 返回行值

}

}

for(j=1;j<=n+1;j++)//交换两行元素{

float s2=a[k][j];

a[k][j]=a[l][j];

a[l][j]=s2;

}

for(j=k;j<=n+1;j++)

//算出第 k 行行元素的值{

float sum1=0;

for(q=1;q

sum1+=a[k][q]*a[q][j];

a[k][j]=a[k][j]-sum1;

}

for(i=k+1;i<=n;i++)

//算出第 k 列列元素的值{

float sum2=0;

for(q=1;q

sum2+=a[i][q]*a[q][k];

a[i][k]=(a[i][k]-sum2)/(a[k][k]);

}

}//第 k 步分解结束for(i=1;i<=n;i++)

{

for(j=1;j<=(n+1);j++)

cout<

cout<<'\n';

}//输出改变后的数组

//输出解 '''''''''''''''''''''''''''''''''''''''''''''''''''''''//

float x[10];

for(i=n-1;i>=1;i--)

{

x[n]=a[n][n+1]/a[n][n];

float sum3=0;for(i=1;i<=n;i++)

for(j=i+1;j<=n;j++){

sum3+=a[i][j]*x[j];cout<<'x'<

x[i]=(a[i][n+1]-sum3)/a[i][i];}//输出解向量

}//回代过程}

结果：方程的解为：

x1= -0.28923 ， x2= 0.34544 ，x3= -0.71281 ，

x4= -0.22061 ， x5= -0.43040 ， x6= 0.15431 ，

x7= -0.057823 ， x8= 0.20105 ，x9= 0.29023 。

分析：我感觉是提高了查错误点的能力和编循环语句的能力，即利用很规整的迭代公式进行

编程。另外列主元三角分解法的阶梯步骤有了更深的了解!

36．逐次超松弛迭代法

x(k)x( k 1)）；

（1）编制解 n 阶线性方程组Ax=b 的 SOR方法的通用程序（要求

（ 2）对于 35 题中所给的线性方程组，取松弛因子i i / 50(i1,2,,99)，容许误差1

10 5，打印松弛因子、迭代次数、最佳松弛因子及解向量。

2

程序为：

#includedo

#include{

#define eps 0.5e-5 // 迭代误差for(i=0;i<9;i++)

void main(void){x0[i]=x1[i];}

{for(i=0;i<9;i++)

int i,j,l;{

float w,t;sum=0;

float m[9];for(j=0;j

float sum;{

float sum=sum+a[i][j]*x1[j];

a[9][9]={{31,-13,0,0,0,-10,0,0,0},{-13,35,-9,0}

,-11,0,0,0,0},{0,-9,31,-10,0,0,0,0,0},{0,0,-10,7for(j=i+1;j<9;j++)

9,-30,0,0,0,-9},{0,0,0,-30,57,-7,0,-5,0},{0,0,0,{sum=sum+a[i][j]*x0[j];}

0,-7,47,-30,0,0},{0,0,0,0,0,-30,41,0,0},{0,0,0,x1[i]=(1-w)*x0[i]+w*(b[i]-sum)/a[i][i];}

0,-5,0,0,27,-2},// 解出九个解

{0,0,0,-9,0,0,0,-2,29}};for(i=0;i<9;i++)

float b[9]={-15,27,-23,0,-20,12,-7,7,10};{

float max(float m[9]);m[i]=x1[i]-x0[i];}

for(t=1;t<=99;t++)l++;

{}while(max(m)>=eps);

l=0;if(max(m)<=eps)

float x0[9]={1,1,1,1,1,1,1,1,1};{cout<<"迭代次数

float x1[9]={1,1,1,1,1,1,1,1,1};="<

w=t/50;for(i=0;i<9;i++)

cout<<"x1"<<"["<

cout<<"-------------------------------------"<<'\n';for(int i=1;i<9;i++)

}{

}if(fabs(m[i])>k)

}k=fabs(m[i]);

float max(float m[9])// 求出最大}

的迭代误差return k;

{}

float k;

结果为：

ω迭代次数ω迭代次数ω迭代次数ω 迭代次数

0.0212910.5255 1.0219 1.5224

0.047000.5453 1.0418 1.5425

0.064860.5651 1.0617 1.5627

0.083730.5848 1.0816 1.5830

0.103030.6046 1.1015 1.6031

0.122550.6244 1.1215 1.6234

0.142210.6442 1.1414 1.6436

0.161940.6641 1.1612 1.6639

0.181730.6839 1.1810 1.6844

0.201560.7037 1.2011 1.7049

0.221420.7236 1.2212 1.7255

0.241300.7434 1.2412 1.7460

0.261200.7633 1.2613 1.7671

0.281110.7832 1.2813 1.7881

0.301030.8030 1.3014 1.8097

0.32960.8229 1.3215 1.82121

0.34900.8428 1.3415 1.84157

0.36850.8627 1.3616 1.86228

0.38800.8826 1.3817 1.88398

0.40750.9025 1.4018 1.901582

0.42710.9224 1.4219 1.925427

0.44680.9423 1.4419 1.942178

0.46640.9622 1.4620 1.961374

0.48610.9821 1.4821 1.981009

0.5058 1.0020 1.5022

从 1.92 到 2.00 均出现不合理的迭代次数，迭代次数偏大。不是合理的迭代值。

当初值为 x=（ 1,1,1,1,1,1,1,1,1,）T时 ,从上表可以看出 ,最佳松弛因子为ω =1.18,迭代次数仅为 10 次, 方程的解为 :

x1= -0.289231 ，x2= 0.345437 ， x3= -0.712811 ，

x4= -0.220608 ，x5= -0.430400 ， x6= 0.154309 ，

x7= -0.057823 ，x8= 0.201054 ， x9= 0.290229 。

数值分析上机题

姓名：戴载星学号： 040139

习题 4

38.（上机） 3 次条插函数

（1）制求第一型 3 次条插函数的通用程序 ;

（2）已知汽曲型点的数据如下：

012345678910 x i

2.51

3.30

4.04 4.70

5.22 5.54 5.78 5.40 5.57 5.70 5.80

y i

端点条件 y0'=0.8 ，y10'=0.2 。用所制程序求的 3 次条插函数S(x), 并打印出

S(i+0.5)(i=0,1,?9)。

解：通用程序：

#include//求出 h[j] 的

void main(void)for(j=0;j<=n-1;j++)

{{

float x[11];//存放数 x[j]h[j]=x[j+1]-x[j];

float y[11];//存放数 y[j]

float h[11];// 存放数 h[j]cout<<'h'<<'['<

float u[11];// 存放数 u[j]}

float v[11];// 存放数 v[j]cout<

float d[11];// 存放数 d[j]//求出 u[j] 和 v[j] 的初

float M[11];//存放数 M[j]v[0]=1;

float b[11];//存放数 b[j]u[n]=1;

float for(j=1;j<=n-1;j++)

t[11],l[11],yy[11],s[4],aa1,aa2,aa3,aa4;{

float s1[10];u[j]=h[j-1]/(h[j-1]+h[j]);

int i,j,n;v[j]=h[j]/(h[j-1]+h[j]);

float xx;//x区}

//将初初始化//求出 d[j] 的

cout<<" 入 n 的 :\n";for(j=1;j

cin>>n;{d[j]=6*((y[j+1]-y[j])/h[j]-(y[j]-y[j-1])/h[

cout<<" 入数 x:\n";j-1])/(h[j]+h[j-1]);}

for(i=0;i<=n;i++)cin>>x[i];d[0]=6*((y[1]-y[0])/h[0]-df[0])/h[0];

cout<<" 入数 y:\n";d[n]=6*(df[1]-(y[n]-y[n-1])/h[n-1])/h[n-1

for(i=0;i<=n;i++)cin>>y[i];];

//入端点for(j=1;j<=n;j++)

float df[2];{

cout<<" 入两个端点 :\n";

for(i=0;i<2;i++)cin>>df[i];cout<<'u'<<'['<

//利用本上的算法求出所需要的}

cout<

for(j=0;j

{}

//将 M[j] 的值输出

cout<<'v'<<'['<

}{cout<<'M'<<'['<

for(j=0;j<=n;j++)//输出插值多项式的系数

{for(j=0;j

{

cout<<'d'<<'['<

cout<

//利用书本上的追赶法求解方程组for(i=0;i<=n;i++)

{b[i]=2;}

cout<

t[0]=b[0];

yy[0]=d[0];

//消元过程

for(i=1;i<=n;i++)

{

l[i]=u[i]/t[i-1];

t[i]=b[i]-l[i]*v[i-1];

yy[i]=d[i]-l[i]*yy[i-1];

}

//回代过程

M[n]=yy[n]/t[n];

for(i=n-1;i>=0;i--)

s[0]=y[j];

s[1]=(y[j+1]-y[j])/h[j]-(M[j]/3+M[j+1]/6)

*h[j];

s[2]=M[j]/2;

s[3]=(M[j+1]-M[j])/(6*h[j]);

cout<<"当x的值在区间"<<'x'<<'['<

的系数 :\n";

for(int k=0;k<4;k++)

{

cout<<'s'<<'['<

}

cout<

}

}

（2）编制的程序求车门的 3 次样条插值函数S(x) ：

x 属于区间 [0,1] 时； S(x)=2.51+0.8(x)-0.0014861(x)(x)-0.00851395(x)(x)(x)

x 属于区间 [1,2] 时； S(x)=3.3+0.771486(x-1)-0.027028(x-1)(x-1)-0.00445799(x-1)(x-1)(x-1)

x 属于区间 [2,3] 时； S(x)=4.04+0.704056(x-2)-0.0404019(x-2)(x-2)-0.0036543(x-2)(x-2)(x-2) x 属于区间 [3,4] 时； S(x)=4.7+0.612289(x-3)-0.0513648(x-3)(x-3)-0.0409245(x-3)(x-3)(x-3)

x 属于区间 [4,5] 时； S(x)=5.22+0.386786(x-4)-0.174138(x-4)(x-4)+0.107352(x-4)(x-4)(x-4)

x 属于区间 [5,6] 时； S(x)=5.54+0.360567(x-5)+0.147919(x-5)(x-5)-0.268485(x-5)(x-5)(x-5)

x 属于区间 [6,7] 时； S(x)=5.78-0.149051(x-6)-0.657537(x-6)(x-6)+0.426588(x-6)(x-6)(x-6)

x 属于区间 [7,8] 时； S(x)=5.4-0.184361(x-7)+0.622227(x-7)(x-7)-0.267865(x-7)(x-7)(x-7)

x 属于区间 [8,9] 时； S(x)=5.57+0.256496(x-8)-0.181369(x-8)(x-8)+0.0548728(x-8)(x-8)(x-8)

x 属于区间 [9,10] 时； S(x)=5.7+0.058376(x-9)-0.0167508(x-9)(x-9)+0.0583752(x-9)(x-9)(x-9) S(0.5)=2.90856S(1.5)=3.67843S (2.5)=4.38147 S(3.5)=4.98819S(4.5)=5.38328S(5.5)=5.7237 S(6.5)=5.59441S(7.5)=5.42989S(8.5)=5.65976 S(9.5)=5.7323

习题五重积分的计算

23（上机题）重积分的计算

d b

题目：给定积分 I ( f )( f (x, y)dx )dy 。取初始步长h和k，及精度。应用复化Simpson

c a

公式，采用逐次二分步长的方法，编制计算I(f) 的通用程序。计算至相邻两次近似值之差的

绝对值不超过为止。

1）用所编程序计算积分I ( f )6 ( 3 tg ( x2y 2 )dx)dy ，取0.5* 10 5。

00

算法概述

初始时候只在 x，y 方向上各二分一次，根据复化 Simpson 公式计算积分值，然后再二分一次，仍然根据上式重新计算积分值，比较两次计算结果的差值，如果小于误差限，则已求得满足要

求的结果，否则继续二分区间直到满足误差要求为止。

程序如下：

#include

#include

#define PI 3.1415926

#define error 0.5E-5

double f(double x,double y)

{

double answ;

answ=tan(x*x+y*y);

return answ;

}

void main()

{

double old,temp,a=0,c=0;

int i,j,m=1,n=1;

double b=PI/3;// 重积分内层上限

double d=PI/6;// 重积分外层上限

double h=(b-a)/(2*n);//

double k=(d-c)/(2*m);

double answ=0;

do

{

old=answ;

answ=0;

answ+=f(a,c);

temp=0;

for(i=1;i

temp+=f(a+2*i*h,c);

temp*=2;

answ+=temp;

temp=0;

for(i=1;i<=n;i++)

temp+=f(a+(2*i-1)*h,c);

temp*=4;

answ+=temp;

answ+=f(b,c);

temp=0;

for(j=1;j

temp+=f(a,c+2*j*k);

temp*=2;

answ+=temp;

temp=0;

for(j=1;j

for(i=1;i

temp+=f(a+2*i*h,c+2*j*k);

temp*=4;

answ+=temp;

temp=0;

for(j=1;j

for(i=1;i<=n;i++)

temp+=f(a+(2*i-1)*h,c+2*j*k);

temp*=8;

answ+=temp;

temp=0;

for(j=1;j

temp+=f(b,c+2*j*k);

temp*=2;

answ+=temp;

temp=0;

for(j=1;j<=m;j++)

temp+=f(a,c+(2*j-1)*k);

temp*=4;

answ+=temp;

temp=0;

for(j=1;j<=m;j++)

for(i=1;i

temp+=f(a+2*i*h,c+(2*j-1)*k);

temp*=8;

answ+=temp;

temp=0;

for(j=1;j<=m;j++)

for(i=1;i<=n;i++)

temp+=f(a+(2*i-1)*h,c+(2*j-1)*k);

temp*=16;

answ+=temp;

temp=0;

for(j=1;j<=m;j++)

temp+=f(b,c+(2*j-1)*k);

temp*=4;

answ+=temp;

answ+=f(a,d);

temp=0;

for(i=1;i

temp+=f(a+2*i*h,d);

temp*=2;

answ+=temp;

temp=0;

for(i=1;i<=n;i++)

temp+=f(a+(2*i-1)*h,d);

temp*=4;

answ+=temp;

answ+=f(b,d);

answ=answ*h*k/9;

m*=2;

n*=2;

h=(b-a)/(2*n);

k=(d-c)/(2*m);

}while(fabs(answ-old)>error);

cout<<"answ is:"<

cout<<"\nDivided into"<

return;

}

程序的输出结果为：

Result is:0.33652

Divided into64parts

即在 x， y 方向上各二分 6 次。本题公式较长，形式较为复杂，但是并不需要太多编程

技巧，只需细心即可。随着二分的继续，计算的结果越来越趋向准确值。本题在二分 6 次后，经检验误差符合精度要求。

习题 6

21．（上机题）常微分方程初值问题数值解

（1）编制 RK 4方法的通用程序；

（2）编制 AB 4方法的通用程序（由RK 4提供初值）；

（3）编制 AB 4 -AM 4预测校正方法的通用程序（由RK4提供初值）；

（4）编制带改进的AB 4-AM 4预测校正方法的通用程序（由RK 4提供初值） ;

（5）对于初值问题

y' x2 y2(0 x 1.5)

y(0) 3

取步长 h0.1 ,应用（1）~（4）中的四种方法进行计算，并将计算结果和精确解y( x) 3/(1x3 ) 作比较；

（6）通过本上机题，你能得到哪些结论？

解：程序为：

#include

#include

#include

#include

ofstream outfile("data.txt");

//此处定义函数f(x,y) 的表达式

//用户可以自己设定所需要求得函数表达式double f1(double x,double y)

{

double f1;

f1=(-1)*x*x*y*y;

return f1;

}

//此处定义求函数精确解的函数表达式double f2(double x)

{

double f2;

f2=3/(1+x*x*x);

return f2;

}

//此处为精确求函数解的通用程序

void accurate(double a,double b,double h){

double x[100],accurate[100];

x[0]=a;

int i=0;

outfile<<" 输出函数准确值的程序结果 :\n";

do{

x[i]=x[0]+i*h;

accurate[i]=f2(x[i]);

outfile<<"accurate["<

i++;

}while(i<(b-a)/h+1);

}

//此处为经典 Runge-Kutta 公式的通用程序void RK4(double a,double b,double h,double c)

{

int i=0;

double k1,k2,k3,k4;

double x[100],y[100];

y[0]=c;

x[0]=a;

outfile<<" 输出经典 Runge-Kutta 公式的

程序结果 :\n";

do y[i+1]=y[i]+h*(k1+2*k2+2*k3+k4)/6;

{}

x[i]=x[0]+i*h;int j=3;

k1=f1(x[i],y[i]);y1[0]=y[0];

k2=f1((x[i]+h/2),(y[i]+h*k1/2));y1[1]=y[1];

k3=f1((x[i]+h/2),(y[i]+h*k2/2));y1[2]=y[2];

k4=f1((x[i]+h),(y[i]+h*k3));y1[3]=y[3];

do

y[i+1]=y[i]+h*(k1+2*k2+2*k3+k4)/6;{

x[j]=x[0]+j*h;

outfile<<"y"<<"["<

i++;y1[j+1]=y1[j]+(55*f1(x[j],y1[j])-59*f1(x }while(i<(b-a)/h+1);[j-1],y1[j-1])+37*f1(x[j-2],y1[j-2])-9*f1(x[j-3 }],y1[j-3]))*h/24;

//此处为 4 阶 Adams 显式方法的通用程序

void AB4(double a,double b,double h,double n';

c)j++;

{}while(j<(b-a)/h+1);

double x[100],y[100],y1[100];}

double k1,k2,k3,k4;

y[0]=c;// 主函数

x[0]=a;void main(void)

outfile<<" 输出 4 阶 Adams 显式方法的{

程序结果 :\n";double a,b,h,c;

for(int i=0;i<=2;i++)cout<<" 输入上下区间、步长和初始{值 :\n";

x[i]=x[0]+i*h;cin>>a>>b>>h>>c;

k1=f1(x[i],y[i]);accurate(a,b,h);

k2=f1((x[i]+h/2),(y[i]+h*k1/2));RK4(a,b,h,c);

k3=f1((x[i]+h/2),(y[i]+h*k2/2));AB4(a,b,h,c);

k4=f1((x[i]+h),(y[i]+h*k3));}

结果为：

由经典 Runge-Kutta 公式得出的结果列在下面的表格中，以及精确值 y(x i )和精确值和数值解的误差：

i x i y i y(x i )|y(x i )-y i |

00330

10.1 2.997 2.997 1.87138e-007

20.2 2.97619 2.97619 3.91665e-007

30.3 2.92113 2.921137.58342e-007

40.4 2.81955 2.81955 1.61101e-006

50.5 2.66666 2.66667 3.17735e-006

60.6 2.4671 2.46711 5.00551e-006

70.7 2.2338 2.2338 5.77233e-006

80.8 1.98412 1.98413 4.12954e-006

90.9 1.73511 1.73511 1.15554e-007

10 1.0 1.50001 1.5 5.80668e-006

11 1.1 1.28701 1.287 1.13075e-005

12 1.2 1.09972 1.09971 1.54242e-005

13 1.30.9383970.93838 1.77272e-005

14 1.40.80130.801282 1.83754e-005

15 1.50.6857320.685714 1.78e-005

由 AB4方法得出的结果为：

Y1[0]=3y1[1]=2.997y1[2]=2.97619y1[3]=2.92113y1[4]=2.81839

y1[5]=2.66467y1[6]=2.4652y1[7]=2.23308y1[8]=1.98495y1[9]=1.73704

y1[10]=1.50219y1[11]=1.28876 y1[12]=1.10072y1[13]=0.93871y1[14]=0.801135

y1[15]=0.685335

通过本上机题我明白了各种求微分方程的数值方法，经典 Runge-Kutta 公式， AB4 方法以及AB4-AM4预测校正方法求解公式的精度是不同的。其中经典Runge-Kutta 公式的精度，四阶 Adams 显式方法（ AB4 ）具有 4 阶精度。

习题八抛物线方程Crank-Nicolson 格式

题目： 1）编制用Crank-Nicolson格式求抛物线方程

u/ t - a2u/2x = f(x,t)(0

u(x,0) =( x)(0x1)

u(0,t) =x，u（1,t）=t(0

数值解的通用程序。

2）就 a=1, f(x,t)=0,x =exp(x),t =exp(t),t =exp(1+t),M=40,N=40,输入点（0.2，

1.0 ），（ 0.4 ， 1.0 ），（ 0.6 ， 1.0 ），（ .8 ， 1.0 ） 4 点处 u（ x,t ）的近似值。

3）已知所给方程的精确解为u(x,t)=exp(x+t),将步长反复二分，从（0.2 ， 1.0 ），（0.4 ，1.0 ），（ 0.6 ， 1.0 ），（ 0.8 ， 1.0 ） 4 点处精确解与数值解的误差观察当步长缩小一半时，

误差以什么规律缩小。

算法：

本题选择空间步长h=0.025, 时间步长ζ =0.025

1.置初值；

2.通过点 (i-1,k),(i-1,k+1),(i,k),(i+1,k),(i+1,k+1)求解点(I,k+1)

的值；

3.输出结果

原程序： #include

#include

float h=0.025,k=0.025;

int m=40;int n=40;

float y[40][40],r=a*k/(h*h);

void Input()

{

int i,j;

cout<<"Loading Input Data..."<

for(i=0;i

{

for(j=0;j

if (i==j) a[i][j]=1+r;

for(j=0;j

if ((j=i+1)||(i=j+1)) a[i][j]=-r/2;

}

}

int main()

{

Input(); //read data

int r,i;

for(k=0;k<(m-1);k++)

{

int select(); //select main element

r=select();

void exchange(int g);

exchange(r); //exchange

void analyze();

analyze(); //analyze

}

void ret();

ret(); // replace back

cout<<"The solution vector is below:"<

cout<<"x["<

return 0;

}

int select()

{

int f,t; float max;

f=k;

float si(int u,int v);

max=float(fabs(si(k,k)));

for(t=(k+1);t<(m-1);t++)

{

if(max

{

max=float(fabs(si(t,k)));

f=t;

}

}

return f;

}

float si(int u,int v)

{

float sum=0; int q;

for(q=0;q

sum+=a[u][q]*a[q][v];

sum=a[u][v]-sum;

return sum;

}

void exchange(int g)

{

int t; float temp;

for(t=0;t

{

temp=a[k][t];

a[k][t]=a[g][t];

a[g][t]=temp;

}

}

void analyze()

{

int t;

float si(int u,int v);

for(t=k;t

a[k][t]=si(k,t);

for(t=(k+1);t

a[t][k]=(float)(si(t,k)/a[k][k]);

}

void ret()

{

int t,z;float sum;

x[m-1]=(float)a[m-1][m]/a[m-1][m-1];

for(t=(m-2);t>-1;t--)

{

sum=0;

for(z=(t+1);z

sum+=a[t][z]*x[z];

x[t]=(float)(a[t][m]-sum)/a[t][t];

}

}

运行结果: ===Numerical Solution of Partial Differential Equations=== NumericalResult AccurateValue Errors

（0.2,1.0 ） 3.31947 3.320120.00045

（0.4,1.0 ） 4.05432 4.055200.00088

（0.6,1.0 ） 4.95238 4.953030.00065

（0.8,1.0 ） 4.04897 6.049650.00067

分析与结论: Crank-Nicolson 格式，步长越小，精度高。但计算较麻烦。

数值分析上机作业

数值分析上机实验报告 选题：曲线拟合的最小二乘法 指导老师： 专业： 学号： 姓名：

课题八曲线拟合的最小二乘法 一、问题提出 从随机的数据中找出其规律性，给出其近似表达式的问题，在生产实践和科学实验中大量存在，通常利用数据的最小二乘法求得拟合曲线。 在某冶炼过程中，根据统计数据的含碳量与时间关系，试求含碳量y 与时间t 的拟合曲线。 二、要求 1、用最小二乘法进行曲线拟合； 2、近似解析表达式为()33221t a t a t a t ++=?； 3、打印出拟合函数()t ?，并打印出()j t ?与()j t y 的误差，12,,2,1 =j ； 4、另外选取一个近似表达式，尝试拟合效果的比较； 5、*绘制出曲线拟合图*。 三、目的和意义 1、掌握曲线拟合的最小二乘法； 2、最小二乘法亦可用于解超定线代数方程组； 3、探索拟合函数的选择与拟合精度间的关系。 四、计算公式 对于给定的测量数据(x i ,f i )(i=1,2,…，n )，设函数分布为 ∑==m j j j x a x y 0)()(? 特别的，取)(x j ?为多项式 j j x x =)(? (j=0, 1,…，m )

则根据最小二乘法原理，可以构造泛函 ∑∑==-=n i m j i j j i m x a f a a a H 1 10))((),,,(? 令 0=??k a H (k=0, 1,…，m ) 则可以得到法方程 ???? ??????? ?=????????????????????????),(),(),(),(),(),(),(),(),(),(),(),(1010101111000100m m m m m m m m f f f a a a ????????????????????? 求该解方程组，则可以得到解m a a a ,,,10 ，因此可得到数据的最小二乘解 ∑=≈m j j j x a x f 0)()(? 曲线拟合：实际工作中，变量间未必都有线性关系，如服药后血药浓度与时间的关系；疾病疗效与疗程长短的关系；毒物剂量与致死率的关系等常呈曲线关系。曲线拟合是指选择适当的曲线类型来拟合观测数据，并用拟合的曲线方程分析两变量间的关系。 五、结构程序设计 在程序结构方面主要是按照顺序结构进行设计，在进行曲线的拟合时，为了进行比较，在程序设计中，直接调用了最小二乘法的拟合函数polyfit ，并且依次调用了plot 、figure 、hold on 函数进行图象的绘制，最后调用了一个绝对值函数abs 用于计算拟合函数与原有数据的误差，进行拟合效果的比较。

数值分析实验报告1

实验一误差分析 实验1.1（病态问题） 实验目的：算法有“优”与“劣”之分，问题也有“好”与“坏”之别。对数值方法的研究而言，所谓坏问题就是问题本身对扰动敏感者，反之属于好问题。通过本实验可获得一个初步体会。 数值分析的大部分研究课题中，如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。病态问题要通过研究和构造特殊的算法来解决，当然一般要付出一些代价（如耗用更多的机器时间、占用更多的存储空间等）。 问题提出：考虑一个高次的代数多项式 显然该多项式的全部根为1,2,…,20共计20个，且每个根都是单重的。现考虑该多项式的一个扰动 其中ε（1.1）和（1.221,,,a a 的输出b ”和“poly ε。 （1（2 （3）写成展 关于α solve 来提高解的精确度，这需要用到将多项式转换为符号多项式的函数poly2sym,函数的具体使用方法可参考Matlab 的帮助。 实验过程： 程序： a=poly(1:20); rr=roots(a); forn=2:21 n form=1:9 ess=10^(-6-m);

ve=zeros(1,21); ve(n)=ess; r=roots(a+ve); -6-m s=max(abs(r-rr)) end end 利用符号函数：（思考题一）a=poly(1:20); y=poly2sym(a); rr=solve(y) n

很容易的得出对一个多次的代数多项式的其中某一项进行很小的扰动，对其多项式的根会有一定的扰动的，所以对于这类病态问题可以借助于MATLAB来进行问题的分析。 学号：06450210 姓名：万轩 实验二插值法

数值分析试题及答案汇总

数值分析试题 一、 填空题（2 0×2′） 1. ?? ????-=? ?????-=32,1223X A 设x =是精确值x *=的近似值，则x 有 2 位 有效数字。 2. 若f (x )=x 7－x 3＋1，则f [20,21,22,23,24,25,26,27]= 1 ， f [20,21,22,23,24,25,26,27,28]= 0 。 3. 设，‖A ‖∞＝___5 ____，‖X ‖∞＝__ 3_____， ‖AX ‖∞≤_15_ __。 4. 非线性方程f (x )=0的迭代函数x =?(x )在有解区间满足 |?’(x )| <1 ，则使用该迭代 函数的迭代解法一定是局部收敛的。 5. 区间[a ,b ]上的三次样条插值函数S (x )在[a ,b ]上具有直到 2 阶的连续导数。 6. 当插值节点为等距分布时，若所求节点靠近首节点，应该选用等距节点下牛顿差商 公式的 前插公式 ，若所求节点靠近尾节点，应该选用等距节点下牛顿差商公式的 后插公式 ；如果要估计结果的舍入误差，应该选用插值公式中的 拉格朗日插值公式 。 7. 拉格朗日插值公式中f (x i )的系数a i (x )的特点是：=∑=n i i x a 0)( 1 ；所以当 系数a i (x )满足 a i (x )>1 ，计算时不会放大f (x i )的误差。 8. 要使 20的近似值的相对误差小于%，至少要取 4 位有效数字。 9. 对任意初始向量X (0)及任意向量g ，线性方程组的迭代公式x (k +1)=Bx (k )+g (k =0,1,…)收 敛于方程组的精确解x *的充分必要条件是 ?(B)<1 。 10. 由下列数据所确定的插值多项式的次数最高是 5 。 11. 牛顿下山法的下山条件为 |f(xn+1)|<|f(xn)| 。 12. 线性方程组的松弛迭代法是通过逐渐减少残差r i (i =0,1,…,n )来实现的，其中的残差 r i ＝ (b i -a i1x 1-a i2x 2-…-a in x n )/a ii ，(i =0,1,…,n )。 13. 在非线性方程f (x )=0使用各种切线法迭代求解时，若在迭代区间存在唯一解，且f (x )

数值分析上机作业

昆明理工大学工科研究生《数值分析》上机实验 学院：材料科学与工程学院 专业：材料物理与化学 学号：2011230024 姓名: 郑录 任课教师：胡杰

P277-E1 1.已知矩阵A= 10787 7565 86109 75910 ?? ?? ?? ?? ?? ??,B= 23456 44567 03678 00289 00010 ?? ?? ?? ?? ?? ?? ?? ?? ,错误！未找到引用源。 = 11/21/31/41/51/6 1/21/31/41/51/61/7 1/31/41/51/61/71/8 1/41/51/61/71/81/9 1/51/61/71/81/91/10 1/61/71/81/91/101/11?????????????????? （1）用MA TLAB函数“eig”求矩阵全部特征值。 （2）用基本QR算法求全部特征值（可用MA TLAB函数“qr”实现矩阵的QR分解）。解：MA TLAB程序如下： 求矩阵A的特征值： clear; A=[10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10]; E=eig(A) 输出结果： 求矩阵B的特征值： clear; B=[2 3 4 5 6;4 4 5 6 7;0 3 6 7 8;0 0 2 8 9;0 0 0 1 0]; E=eig(B) 输出结果：

求矩阵错误！未找到引用源。的特征值： clear; 错误！未找到引用源。=[1 1/2 1/3 1/4 1/5 1/6; 1/2 1/3 1/4 1/5 1/6 1/7; 1/3 1/4 1/5 1/6 1/7 1/8; 1/4 1/5 1/6 1/7 1/8 1/9;1/5 1/6 1/7 1/8 1/9 1/10; 1/6 1/7 1/8 1/9 1/10 1/11]; E=eig(错误！未找到引用源。) 输出结果： （2）A= 10 7877565861097 5 9 10 第一步：A0=hess(A);[Q0,R0]=qr(A0);A1=R0*Q0 返回得到： 第二部：[Q1,R1]=qr(A1);A2=R1*Q1

数值分析课后题答案

数值分析 第二章 2．当1,1,2x =-时，()0,3,4f x =-,求()f x 的二次插值多项式。 解： 0120121200102021101201220211,1,2, ()0,()3,()4;()()1 ()(1)(2)()()2()()1 ()(1)(2) ()()6 ()()1 ()(1)(1) ()()3 x x x f x f x f x x x x x l x x x x x x x x x x x l x x x x x x x x x x x l x x x x x x x ==-===-=--==-+-----==------= =-+-- 则二次拉格朗日插值多项式为 2 20 ()()k k k L x y l x ==∑ 0223()4() 14 (1)(2)(1)(1)23 537623 l x l x x x x x x x =-+=---+ -+= +- 6．设,0,1,,j x j n =L 为互异节点，求证： （1） 0()n k k j j j x l x x =≡∑ (0,1,,);k n =L （2） ()()0n k j j j x x l x =-≡∑ (0,1,,);k n =L 证明 （1） 令()k f x x = 若插值节点为,0,1,,j x j n =L ，则函数()f x 的n 次插值多项式为0 ()()n k n j j j L x x l x == ∑。

插值余项为(1)1() ()()()()(1)! n n n n f R x f x L x x n ξω++=-= + 又,k n ≤Q (1)()0 ()0 n n f R x ξ+∴=∴= 0()n k k j j j x l x x =∴=∑ (0,1,,);k n =L 0 000 (2)()() (())()()(()) n k j j j n n j i k i k j j j i n n i k i i k j j i j x x l x C x x l x C x x l x =-==-==-=-=-∑∑∑∑∑ 0i n ≤≤Q 又 由上题结论可知 ()n k i j j j x l x x ==∑ ()()0 n i k i i k i k C x x x x -=∴=-=-=∑原式 ∴得证。 7设[]2 (),f x C a b ∈且()()0,f a f b ==求证： 21 max ()()max ().8 a x b a x b f x b a f x ≤≤≤≤''≤- 解：令01,x a x b ==，以此为插值节点，则线性插值多项式为 10 101010 ()() ()x x x x L x f x f x x x x x --=+-- =() ()x b x a f a f b a b x a --=+-- 1()()0()0 f a f b L x ==∴=Q 又

数值分析上机题课后作业全部-东南大学

2015.1.9 上机作业题报告 USER

1.Chapter1 1.1题目 设S N = 1 j 2?1 N j =2 ，其精确值为 )1 1123(21+--N N 。 （1）编制按从大到小的顺序1 1 1311212 22-+??+-+-=N S N ，计算S N 的通用程序。 （2）编制按从小到大的顺序1 21 1)1(111222-+ ??+--+-= N N S N ，计算S N 的通用程序。 （3）按两种顺序分别计算64210,10,10S S S ,并指出有效位数。（编制程序时用单精度） （4）通过本次上机题，你明白了什么？ 1.2程序 1.3运行结果

1.4结果分析 按从大到小的顺序，有效位数分别为：6，4，3。 按从小到大的顺序，有效位数分别为：5，6，6。 可以看出，不同的算法造成的误差限是不同的，好的算法可以让结果更加精确。当采用从大到小的顺序累加的算法时，误差限随着N 的增大而增大，可见在累加的过程中，误差在放大，造成结果的误差较大。因此，采取从小到大的顺序累加得到的结果更加精确。 2.Chapter2 2.1题目 （1）给定初值0x 及容许误差ε，编制牛顿法解方程f(x)=0的通用程序。 （2）给定方程03 )(3 =-=x x x f ,易知其有三个根3,0,3321= *=*-=*x x x ○1由牛顿方法的局部收敛性可知存在,0>δ当),(0δδ+-∈x 时，Newton 迭代序列收敛于根x2*。试确定尽可能大的δ。 ○2试取若干初始值，观察当),1(),1,(),,(),,1(),1,(0+∞+-----∞∈δδδδx 时Newton 序列的收敛性以及收敛于哪一个根。 （3）通过本上机题，你明白了什么？ 2.2程序

东南大学数值分析上机作业汇总

东南大学数值分析上机作业 汇总 -标准化文件发布号：（9456-EUATWK-MWUB-WUNN-INNUL-DDQTY-KII

数值分析上机报告 院系： 学号： 姓名：

目录 作业1、舍入误差与有效数 (1) 1、函数文件cxdd.m (1) 2、函数文件cddx.m (1) 3、两种方法有效位数对比 (1) 4、心得 (2) 作业2、Newton迭代法 (2) 1、通用程序函数文件 (3) 2、局部收敛性 (4) （1）最大δ值文件 (4) （2）验证局部收敛性 (4) 3、心得 (6) 作业3、列主元素Gauss消去法 (7) 1、列主元Gauss消去法的通用程序 (7) 2、解题中线性方程组 (7) 3、心得 (9) 作业4、三次样条插值函数 (10) 1、第一型三次样条插值函数通用程序： (10) 2、数据输入及计算结果 (12)

作业1、舍入误差与有效数 设∑ =-=N j N j S 2 2 11 ，其精确值为?? ? ??---1112321N N . （1）编制按从小到大的顺序1 1 131121222-? ??+-+-=N S N ，计算N S 的通用程序； （2）编制按从大到小的顺序()1 21 11111222-???+--+-=N N S N ，计算N S 的通用程序； （3）按两种顺序分别计算642101010,,S S S ,并指出有效位数； （4）通过本上机你明白了什么？ 程序： 1、函数文件cxdd.m function S=cxdd(N) S=0; i=2.0; while (i<=N) S=S+1.0/(i*i-1); i=i+1; end script 运行结果（省略>>）： S=cxdd(80) S= 0.737577 2、函数文件cddx.m function S=cddx (N) S=0; for i=N:-1:2 S=S+1/(i*i-1); end script 运行结果（省略>>）： S=cddx(80) S= 0.737577 3、两种方法有效位数对比

数值分析上机实验报告

数值分析上机实验报告

《数值分析》上机实验报告 1.用Newton 法求方程 X 7-X 4+14=0 在（0.1，1.9）中的近似根（初始近似值取为区间端点，迭代6次或误差小于0.00001）。 1.1 理论依据： 设函数在有限区间[a ，b]上二阶导数存在，且满足条件 {}α?上的惟一解在区间平方收敛于方程所生的迭代序列 迭代过程由则对任意初始近似值达到的一个中使是其中上不变号 在区间],[0)(3,2,1,0,) (') ()(],,[x |))(),((|,|,)(||)(|.4;0)(.3],[)(.20 )()(.110......b a x f x k x f x f x x x Newton b a b f a f mir b a c x f a b c f x f b a x f b f x f k k k k k k ==- ==∈≤-≠>+ 令 )9.1()9.1(0)8(4233642)(0)16(71127)(0)9.1(,0)1.0(,1428)(3 2 2 5 333647>?''<-=-=''<-=-='<>+-=f f x x x x x f x x x x x f f f x x x f 故以1.9为起点 ?? ?? ? ='- =+9.1)()(01x x f x f x x k k k k 如此一次一次的迭代，逼近x 的真实根。当前后两个的差<=ε时，就认为求出了近似的根。本程序用Newton 法求代数方程(最高次数不大于10)在（a,b ）区间的根。

1.2 C语言程序原代码： #include #include main() {double x2,f,f1; double x1=1.9; //取初值为1.9 do {x2=x1; f=pow(x2,7)-28*pow(x2,4)+14; f1=7*pow(x2,6)-4*28*pow(x2,3); x1=x2-f/f1;} while(fabs(x1-x2)>=0.00001||x1<0.1); //限制循环次数printf("计算结果：x=%f\n",x1);} 1.3 运行结果： 1.4 MATLAB上机程序 function y=Newton(f,df,x0,eps,M) d=0; for k=1:M if feval(df,x0)==0 d=2;break else x1=x0-feval(f,x0)/feval(df,x0); end e=abs(x1-x0); x0=x1; if e<=eps&&abs(feval(f,x1))<=eps d=1;break end end

数值分析试题及答案

一、单项选择题（每小题3分，共15分） 1. 3.142和3.141分别作为π的近似数具有（ ）和（ ）位有效数字. A ．4和3 B ．3和2 C ．3和4 D ．4和4 2. 已知求积公式 ()()2 1 121 1()(2)636f x dx f Af f ≈ ++? ，则A ＝（ ） A ． 16 B ．13 C ．12 D ．2 3 3. 通过点 ()()0011,,,x y x y 的拉格朗日插值基函数()()01,l x l x 满足（ ） A ． ()00l x ＝0， ()110l x = B ． ()00l x ＝0， ()111l x = C ．() 00l x ＝1，()111 l x = D ． () 00l x ＝1，()111 l x = 4. 设求方程 ()0 f x =的根的牛顿法收敛，则它具有（ ）敛速。 A ．超线性 B ．平方 C ．线性 D ．三次 5. 用列主元消元法解线性方程组 1231231 220223332 x x x x x x x x ++=?? ++=??--=? 作第一次消元后得到的第3个方程（ ）. A ． 232 x x -+= B ．232 1.5 3.5 x x -+= C ． 2323 x x -+= D ． 230.5 1.5 x x -=- 单项选择题答案 1.A 2.D 3.D 4.C 5.B 得 分 评卷人 二、填空题（每小题3分，共15分）

1. 设T X )4,3,2(-=, 则=1||||X ，2||||X = . 2. 一阶均差 ()01,f x x = 3. 已知3n =时，科茨系数()()() 33301213,88C C C ===，那么 () 33C = 4. 因为方程()420 x f x x =-+=在区间 []1,2上满足 ，所以()0f x =在区间 内有根。 5. 取步长0.1h =，用欧拉法解初值问题 ()211y y y x y ?'=+?? ?=? 的计算公式 . 填空题答案 1. 9和29 2. ()() 0101 f x f x x x -- 3. 1 8 4. ()()120 f f < 5. ()12 00.1 1.1,0,1,210.11k k y y k k y +???? ?=+? ?=+???? =??L 得 分 评卷人 三、计算题（每题15分，共60分） 1. 已知函数 21 1y x = +的一组数据： 求分 段线性插值函数，并计算 () 1.5f 的近似值. 计算题1.答案 1. 解 []0,1x ∈， ()1010.510.50110x x L x x --=?+?=---% []1,2x ∈，()210.50.20.30.81221x x L x x --=?+?=-+--%

数值分析作业思考题汇总

￥ 数值分析思考题1 1、讨论绝对误差（限）、相对误差（限）与有效数字之间的关系。 2、相对误差在什么情况下可以用下式代替 3、查阅何谓问题的“病态性”，并区分与“数值稳定性”的不同点。 4、取 ，计算 ，下列方法中哪种最好为什么（1）(3 3-，（2）(2 7-，（3） ()3 1 3+ ，（4） ()6 1 1 ，（5）99- , 数值实验 数值实验综述：线性代数方程组的解法是一切科学计算的基础与核心问题。求解方法大致可分为直接法和迭代法两大类。直接法——指在没有舍入误差的情况下经过有限次运算可求得方程组的精确解的方法，因此也称为精确法。当系数矩阵是方的、稠密的、无任何特殊结构的中小规模线性方程组时，Gauss消去法是目前最基本和常用的方法。如若系数矩阵具有某种特殊形式，则为了尽可能地减少计算量与存储量，需采用其他专门的方法来求解。 Gauss消去等同于矩阵的三角分解，但它存在潜在的不稳定性，故需要选主元素。对正定对称矩阵，采用平方根方法无需选主元。方程组的性态与方程组的条件数有关，对于病态的方程组必须采用特殊的方法进行求解。 数值计算方法上机题目1 1、实验1. 病态问题 实验目的： 算法有“优”与“劣”之分，问题也有“好”和“坏”之别。所谓坏问题就是问题本身的解对数据变化的比较敏感，反之属于好问题。希望读者通过本实验对此有一个初步的体会。 数值分析的大部分研究课题中，如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。病态问题要通过研究和构造特殊的算法来解决，当然一般要付出一些代价（如耗用更多的机器时间、占用更多的存储空间等）。 $ r e x x e x x ** * ** - == 141 . ≈)61

Matlab作业3(数值分析)答案

Matlab作业3（数值分析） 机电工程学院（院、系）专业班组 学号姓名实验日期教师评定 1.计算多项式乘法(x2+2x+2)(x2+5x+4)。 答： 2. （1）将(x-6)(x-3)(x-8)展开为系数多项式的形式。（2）求解在x=8时多项 式(x-1)(x-2) (x-3)(x-4)的值。 答：（1） （2）

3. y=sin(x)，x从0到2π，?x=0.02π，求y的最大值、最小值、均值和标准差。 4.设ｘ＝[0.00.30.8 1.1 1.6 2.3]'，y=[0.500.82 1.14 1.25 1.35 1.40]'，试求二次多项式拟合系数，并据此计算x1=[0.9 1.2]时对应的y1。解：x=[0.0 0.3 0.8 1.1 1.6 2.3]'; %输入变量数据x y=[0.50 0.82 1.14 1.25 1.35 1.40]'; %输入变量数据y p=polyfit(x,y,2) %对x,y用二次多项式拟合,得到系数p x1=[0.9 1.2]; %输入点x1 y1=polyval(p,x1) %估计x1处对应的y1 p = -0.2387 0.9191 0.5318 y1 = a) 1.2909

5.实验数据处理：已知某压力传感器的测试数据如下表 p为压力值，u为电压值，试用多项式 d cp bp ap p u+ + + =2 3 ) ( 来拟 合其特性函数，求出a,b,c,d，并把拟合曲线和各个测试数据点画在同一幅图上。解： >> p=[0.0,1.1,2.1,2.8,4.2,5.0,6.1,6.9,8.1,9.0,9.9]; u=[10,11,13,14,17,18,22,24,29,34,39]; x=polyfit(p,u,3) %得多项式系数 t=linspace(0,10,100); y=polyval(x,t); %求多项式得值 plot(p,u,'*',t,y,'r') %画拟和曲线 x = 0.0195 -0.0412 1.4469 9.8267

数值分析拉格朗日插值法上机实验报告

课题一：拉格朗日插值法 1.实验目的 1．学习和掌握拉格朗日插值多项式。 2.运用拉格朗日插值多项式进行计算。 2.实验过程 作出插值点（1.00，0.00），（-1.00，-3.00），（2.00，4.00）二、算法步骤 已知：某些点的坐标以及点数。 输入：条件点数以及这些点的坐标。 输出：根据给定的点求出其对应的拉格朗日插值多项式的值。 3.程序流程： （1）输入已知点的个数； （2）分别输入已知点的X坐标； （3）分别输入已知点的Y坐标； 程序如下： #include #include #include float lagrange(float *x,float *y,float xx,int n) /*拉格朗日

插值算法*/ { int i,j; float *a,yy=0.0; /*a作为临时变量，记录拉格朗日插值多项*/ a=(float*)malloc(n*sizeof(float)); for(i=0;i<=n-1;i++) { a[i]=y[i]; for(j=0;j<=n-1;j++) if(j!=i) a[i]*=(xx-x[j])/(x[i]-x[j]); yy+=a[i]; } free(a); return yy; } int main() { int i; int n; float x[20],y[20],xx,yy; printf("Input n:");

scanf("%d",&n); if(n<=0) { printf("Error! The value of n must in (0,20)."); getch();return 1; } for(i=0;i<=n-1;i++) { printf("x[%d]:",i); scanf("%f",&x[i]); } printf("\n"); for(i=0;i<=n-1;i++) { printf("y[%d]:",i);scanf("%f",&y[i]); } printf("\n"); printf("Input xx:"); scanf("%f",&xx); yy=lagrange(x,y,xx,n); printf("x=%f,y=%f\n",xx,yy); getch(); } 举例如下：已知当x=1，-1,2时f(x)=0，-3,4，求f(1.5)的值。

数值分析整理版试题及答案

数值分析整理版试题及答案

例1、 已知函数表 x -1 1 2 ()f x -3 0 4 求()f x 的Lagrange 二次插值多项式和Newton 二次插值多项式。 解： （1）k x -1 1 2 k y -3 0 4 插值基函数分别为 ()()()()()()()()()() 1200102121()1211126 x x x x x x l x x x x x x x ----= ==-------- ()()()()()()()() ()()021******* ()1211122x x x x x x l x x x x x x x --+-= ==-+---+- ()()()()()()()()()()0122021111 ()1121213 x x x x x x l x x x x x x x --+-= ==-+--+- 故所求二次拉格朗日插值多项式为 () ()()()()()()()()()()2 20 2()11131201241162314 121123537623k k k L x y l x x x x x x x x x x x x x ==?? =-? --+?-+-+?+-????=---++-=+-∑ （2）一阶均差、二阶均差分别为

[]()()[]()()[][][]010********* 011201202303 ,11204 ,412 3 4,,5 2,,126 f x f x f x x x x f x f x f x x x x f x x f x x f x x x x x ---===-----= = =----=== --- k x ()k f x 一阶 二阶 -1 -3 1 0 3/ 2 2 4 4 5/6 故所求Newton 二次插值多项式为 ()()[]()[]()() ()()()20010012012,,,35 311126537623P x f x f x x x x f x x x x x x x x x x x x =+-+--=-+ +++-=+- 例2、 设2 ()32f x x x =++，[0,1]x ∈，试求()f x 在[0, 1]上关于()1x ρ=，{} span 1,x Φ=的最佳平方逼近多项式。 解： 若{}span 1,x Φ=，则0()1x ?=，1()x x ?=，且()1x ρ=，这样，有

(完整版)数值计算方法上机实习题答案

1． 设?+=1 05dx x x I n n ， （1） 由递推公式n I I n n 1 51+-=-，从0I 的几个近似值出发，计算20I ； 解：易得：0I =ln6-ln5=0.1823, 程序为： I=0.182; for n=1:20 I=(-5)*I+1/n; end I 输出结果为：20I = -3.0666e+010 （2） 粗糙估计20I ，用n I I n n 51 5111+- =--，计算0I ； 因为 0095.05 6 0079.01020 201 020 ≈<<≈??dx x I dx x 所以取0087.0)0095.00079.0(2 1 20=+= I 程序为：I=0.0087; for n=1:20 I=(-1/5)*I+1/(5*n); end I 0I = 0.0083 （3） 分析结果的可靠性及产生此现象的原因(重点分析原因)。 首先分析两种递推式的误差；设第一递推式中开始时的误差为000I I E '-=，递推过程的舍入误差不计。并记n n n I I E '-=，则有01)5(5E E E n n n -==-=-Λ。因为=20E 20020)5(I E >>-，所此递推式不可靠。而在第二种递推式中n n E E E )5 1(5110-==-=Λ，误差在缩小， 所以此递推式是可靠的。出现以上运行结果的主要原因是在构造递推式过程中，考虑误差是否得到控制， 即算法是否数值稳定。 2． 求方程0210=-+x e x 的近似根，要求4 1105-+?<-k k x x ，并比较计算量。 （1） 在[0，1]上用二分法； 程序：a=0;b=1.0; while abs(b-a)>5*1e-4 c=(b+a)/2;

数值分析实验报告1

实验一 误差分析 实验（病态问题） 实验目的：算法有“优”与“劣”之分，问题也有“好”与“坏”之别。对数值方法的研究而言，所谓坏问题就是问题本身对扰动敏感者，反之属于好问题。通过本实验可获得一个初步体会。 数值分析的大部分研究课题中，如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。病态问题要通过研究和构造特殊的算法来解决，当然一般要付出一些代价（如耗用更多的机器时间、占用更多的存储空间等）。 问题提出：考虑一个高次的代数多项式 )1.1() ()20()2)(1()(20 1∏=-=---=k k x x x x x p 显然该多项式的全部根为1,2,…,20共计20个，且每个根都是单重的。现考虑该多项式的一个扰动 )2.1(0 )(19=+x x p ε 其中ε是一个非常小的数。这相当于是对（）中19x 的系数作一个小的扰动。我们希望比较（）和（）根的差别，从而分析方程（）的解对扰动的敏感性。 实验内容：为了实现方便，我们先介绍两个Matlab 函数：“roots ”和“poly ”。 roots(a)u = 其中若变量a 存储n+1维的向量，则该函数的输出u 为一个n 维的向量。设a 的元素依次为121,,,+n a a a ，则输出u 的各分量是多项式方程 01121=+++++-n n n n a x a x a x a 的全部根；而函数 poly(v)b =

的输出b 是一个n+1维变量，它是以n 维变量v 的各分量为根的多项式的系数。可见“roots ”和“poly ”是两个互逆的运算函数。 ;000000001.0=ess );21,1(zeros ve = ;)2(ess ve = ))20:1((ve poly roots + 上述简单的Matlab 程序便得到（）的全部根，程序中的“ess ”即是（）中的ε。 实验要求： （1）选择充分小的ess ，反复进行上述实验，记录结果的变化并分析它们。 如果扰动项的系数ε很小，我们自然感觉（）和（）的解应当相差很小。计算中你有什么出乎意料的发现表明有些解关于如此的扰动敏感性如何 （2）将方程（）中的扰动项改成18x ε或其它形式，实验中又有怎样的现象 出现 （3）（选作部分）请从理论上分析产生这一问题的根源。注意我们可以将 方程（）写成展开的形式， ) 3.1(0 ),(1920=+-= x x x p αα 同时将方程的解x 看成是系数α的函数，考察方程的某个解关于α的扰动是否敏感，与研究它关于α的导数的大小有何关系为什么你发现了什么现象，哪些根关于α的变化更敏感 思考题一：（上述实验的改进） 在上述实验中我们会发现用roots 函数求解多项式方程的精度不高，为此你可以考虑用符号函数solve 来提高解的精确度，这需要用到将多项式转换为符号多项式的函数poly2sym,函数的具体使用方法可参考Matlab 的帮助。

数值分析第四版习题及答案

第四版 数值分析习题 第一章 绪 论 1. 设x >0,x 的相对误差为δ,求ln x 的误差. 2. 设x 的相对误差为2％,求n x 的相对误差. 3. 下列各数都是经过四舍五入得到的近似数,即误差限不超过最后一位的半个单位,试指 出它们是几位有效数字: *****123451.1021,0.031,385.6,56.430,7 1.0.x x x x x =====? 4. 利用公式求下列各近似值的误差限: ********12412324(),(),()/,i x x x ii x x x iii x x ++其中**** 1234 ,,,x x x x 均为第3题所给的数. 5. 计算球体积要使相对误差限为1％,问度量半径R 时允许的相对误差限是多少? 6. 设028,Y =按递推公式 1n n Y Y -=…) 计算到100Y .(五位有效数字),试问计算100Y 将有多大误差? 7. 求方程2 5610x x -+=的两个根,使它至少具有四位有效数字. 8. 当N 充分大时,怎样求 2 11N dx x +∞ +? ? 9. 正方形的边长大约为100㎝,应怎样测量才能使其面积误差不超过1㎝2 ? 10. 设 212S gt = 假定g 是准确的,而对t 的测量有±秒的误差,证明当t 增加时S 的绝对误 差增加,而相对误差却减小. 11. 序列 {}n y 满足递推关系1101n n y y -=-(n=1,2,…),若0 1.41y =≈(三位有效数字), 计算到 10y 时误差有多大?这个计算过程稳定吗? 12. 计算61)f =, 1.4≈,利用下列等式计算,哪一个得到的结果最好? 3 -- 13. ()ln(f x x =,求f (30)的值.若开平方用六位函数表,问求对数时误差有多大?若改用另一等价公式 ln(ln(x x =- 计算,求对数时误差有多大?

数值分析实验报告模板

数值分析实验报告模板 篇一：数值分析实验报告(一)(完整) 数值分析实验报告 1 2 3 4 5 篇二：数值分析实验报告 实验报告一 题目：非线性方程求解 摘要：非线性方程的解析解通常很难给出，因此线性方程的数值解法就尤为重要。本实验采用两种常见的求解方法二分法和Newton法及改进的Newton法。利用二分法求解给定非线性方程的根，在给定的范围内，假设f(x,y)在[a,b]上连续，f(a)xf(b) 直接影响迭代的次数甚至迭代的收敛与发散。即若x0 偏离所求根较远，Newton法可能发散的结论。并且本实验中还利用利用改进的Newton法求解同样的方程，且将结果与Newton法的结果比较分析。 前言：（目的和意义） 掌握二分法与Newton法的基本原理和应用。掌握二分法的原理，验证二分法，在选对有根区间的前提下，必是收

敛，但精度不够。熟悉Matlab语言编程，学习编程要点。体会Newton使用时的优点，和局部收敛性，而在初值选取不当时，会发散。 数学原理： 对于一个非线性方程的数值解法很多。在此介绍两种最常见的方法：二分法和Newton法。 对于二分法，其数学实质就是说对于给定的待求解的方程f(x)，其在[a,b]上连续，f(a)f(b) Newton法通常预先要给出一个猜测初值x0，然后根据其迭代公式xk?1?xk?f(xk) f'(xk) 产生逼近解x*的迭代数列{xk}，这就是Newton法的思想。当x0接近x*时收敛很快，但是当x0选择不好时，可能会发散，因此初值的选取很重要。另外，若将该迭代公式改进为 xk?1?xk?rf(xk) 'f(xk) 其中r为要求的方程的根的重数，这就是改进的Newton 法，当求解已知重数的方程的根时，在同种条件下其收敛速度要比Newton法快的多。 程序设计： 本实验采用Matlab的M文件编写。其中待求解的方程写成function的方式，如下 function y=f(x);

数值分析习题集及答案Word版

数值分析习题集 （适合课程《数值方法A 》和《数值方法B 》） 长沙理工大学 第一章 绪 论 1. 设x >0,x 的相对误差为δ,求ln x 的误差. 2. 设x 的相对误差为2％,求n x 的相对误差. 3. 下列各数都是经过四舍五入得到的近似数,即误差限不超过最后一位的半个单位,试指 出它们是几位有效数字: *****123451.1021,0.031,385.6,56.430,7 1.0.x x x x x =====? 4. 利用公式(3.3)求下列各近似值的误差限: ********12412324(),(),()/,i x x x ii x x x iii x x ++其中**** 1234 ,,,x x x x 均为第3题所给的数. 5. 计算球体积要使相对误差限为1％,问度量半径R 时允许的相对误差限是多少? 6. 设028,Y =按递推公式 1n n Y Y -=…) 计算到100Y .27.982(五位有效数字),试问计算100Y 将有多大误差? 7. 求方程2 5610x x -+=的两个根,使它至少具有四位有效数字27.982). 8. 当N 充分大时,怎样求2 1 1N dx x +∞+?? 9. 正方形的边长大约为100㎝,应怎样测量才能使其面积误差不超过1㎝2 ? 10. 设 212S gt = 假定g 是准确的,而对t 的测量有±0.1秒的误差,证明当t 增加时S 的绝对 误差增加,而相对误差却减小. 11. 序列 {}n y 满足递推关系1101n n y y -=-(n=1,2,…),若0 1.41y =≈(三位有效数字), 计算到 10y 时误差有多大?这个计算过程稳定吗? 12. 计算6 1)f =, 1.4≈,利用下列等式计算,哪一个得到的结果最好? 3 -- 13. ()ln(f x x =,求f (30)的值.若开平方用六位函数表,问求对数时误差有多大?

数值分析上机作业1-1

数值计算方法上机题目1 1、实验1. 病态问题 实验目的： 算法有“优”与“劣”之分，问题也有“好”和“坏”之别。所谓坏问题就是问题本身的解对数据变化的比较敏感，反之属于好问题。希望读者通过本实验对此有一个初步的体会。 数值分析的大部分研究课题中，如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。病态问题要通过研究和构造特殊的算法来解决，当然一般要付出一些代价（如耗用更多的机器时间、占用更多的存储空间等）。 问题提出： 考虑一个高次的代数多项式 ∏=-= ---=20 1)()20)...(2)(1()(k k x x x x x p （E1-1） 显然该多项式的全部根为l ，2，…，20，共计20个，且每个根都是单重的（也称为简 单的）。现考虑该多项式方程的一个扰动 0)(19 =+x x p ε (E1-2) 其中ε是一个非常小的数。这相当于是对（E1-1）中19 x 的系数作一个小的扰动。我们希望比较（E1-1）和（E1-2）根的差别，从而分析方程（E1-1）的解对扰动的敏感性。 实验内容： 为了实现方便，我们先介绍两个 Matlab 函数：“roots ”和“poly ”，输入函数 u ＝roots （a ） 其中若变量a 存储1+n 维的向量，则该函数的输出u 为一个n 维的向量。设a 的元素依次为121,...,,+n a a a ，则输出u 的各分量是多项式方程 0...1121=++++-n n n n a x a x a x a 的全部根，而函数 b=poly(v) 的输出b 是一个n ＋1维变量，它是以n 维变量v 的各分量为根的多项式的系数。可见“roots ”和“Poly ”是两个互逆的运算函数. ve=zeros(1,21); ve(2)=ess; roots(poly(1:20))+ve) 上述简单的Matlab 程序便得到（E1-2）的全部根，程序中的“ess ”即是（E1-2）中的ε。 实验要求： （1）选择充分小的ess ，反复进行上述实验，记录结果的变化并分析它们。如果扰动项的系数ε很小，我们自然感觉（E1-1）和(E1-2)的解应当相差很小。计算中你有什么出乎意料的发现？表明有些解关于如此的扰动敏感性如何？ （2）将方程（E1-2）中的扰动项改成18 x ε或其他形式，实验中又有怎样的现象出现？