数值分析实验报告Word文件下载.docx
《数值分析实验报告Word文件下载.docx》由会员分享，可在线阅读，更多相关《数值分析实验报告Word文件下载.docx（13页珍藏版）》请在冰豆网上搜索。
//n整型变量，给定结点的个数
intnewton（doublex[],doubley[],intn）
{
if（n<
=1）/*结点个数不对，此方法不适用，返回*/
return（0）;
else{
for（inti=1;
i<
n;
i++）
{
for（intj=n-1;
j>
=i;
j--）
{
y[j]=（y[j]-y[j-1]）/（x[j]-x[j-i]）;
}
}
return1;
}
}
voidmain（）
{
doublez,s,t=0.2;
//t为插值结点
staticdoublex[5]={0.2,0.4,0.6,0.8,1.0};
staticdoubley[5]={0.98,0.92,0.81,0.64,0.38};
inti,flag=1,n=5;
printf（"
请输入插值点:
"
）;
scanf（"
%lf"
&
t）;
//获取待插值结点
flag=newton（x,y,n）;
//调用求差商的函数，求得牛顿插值多项式中的各阶差商
if（flag==0）
{
printf（"
牛顿插值法不适用于此题！
\n"
else
{//求t点处的牛顿插值
z=y[0];
s=1;
for（i=1;
s=s*（t-x[i-1]）;
//差商前的系数
z+=y[i]*s;
x=%g,y=%g\n"
t,z）;
//输出结果
3、实验结果：
实验二数值积分
1、实验内容
计算积分
，用自适应辛普森积分，使其精度达到0.000001。
分析辛普森算法，其实现如下，其中simp（）实现辛普森算法，fexpression（）为被积函数表达式，main（）为主函数程序，为了测试。
#include<
#include"
math.h"
自适应辛普生积分，函数返回一个双精度实型积分值。
//a双精度实型变量,积分下限。
//b双精度实型变量,积分上限,要求b>
a。
//eps双精度实型变量,积分精度要求。
//f双精度函数指针变量,指向计算被积函数值的函数。
doublesimp（doublea,doubleb,doubleeps,double（*f）（double））
intn,k;
doubleh,t1,t2,s1,s2,ep,p,x;
n=1;
h=b-a;
t1=h*（（*f）（a）+（*f）（b））/2.0;
s1=t1;
/*用T1代替S0*/
ep=eps+1.0;
while（ep>
=eps）
{
p=0.0;
for（k=0;
k<
=n-1;
k++）
x=a+（k+0.5）*h;
p=p+（*f）（x）;
t2=（t1+h*p）/2.0;
/*计算T1=（b-a）*（f（a）+f（b））/2*/
s2=（4.0*t2-t1）/3.0;
/*计算Sn=（4T2n-Tn）/3*/
ep=fabs（s2-s1）;
/*计算精度*/
t1=t2;
s1=s2;
n=n+n;
h=h/2.0;
}
return（s2）;
doublefexpression（doublex）//被积函数表达式y=f（x）
doubley;
y=log（1.0+x）/（1.0+x*x）;
return（y）;
voidmain（）//主函数程序
doublea,b,eps,t,fexpression（double）;
a=0.0;
b=1.0;
eps=0.000001;
t=simp（a,b,eps,fexpression）;
t=%g\n"
t）;
实验三解线性方程组的直接接法
1、实验内容
用LU分解及列主元高斯消去法解线性方程组
=
2、实验分析与实现
（1）LU分解
分析LU分解法解方程组算法，L为对角线元素为1的下三角矩阵，U为上三角矩阵，可用一个矩阵表示。
三角分解法解方程组
//参数说明
//n整型变量。
方程组的阶数
//a双精度实型一维数组，体积为n*n，存放方程组的系数矩阵
//u双精度实型一维数组，体积为n*n,返回时上三角存放U中上三角元素，对角线以下存放L的对角线以下元素
//x双精度实型一维数组，长度为n,存放方程组的中间解和最后的解
//b双精度实型一维数组，长度为n,存放方程组右端的常数向量
voidtriange（doublea[],doubleu[],doublex[],doubleb[],intn）
inti,k,r;
for（i=0;
i++）{
u[0*n+i]=a[0*n+i];
//计算U的第一行元素
for（i=1;
u[i*n+0]=a[i*n+0]/u[0];
//计算L的第一列元素
for（r=1;
r<
r++）{/*计算U中第r行元素，L中第r列元素*/
for（i=r;
i++）{//计算U中元素u[r][i]
u[r*n+i]=a[r*n+i];
for（k=0;
k<
=r-1;
k++）
u[r*n+i]-=u[r*n+k]*u[k*n+i];
for（i=r+1;
n&
&
r!
i++）{//计算L中元素u[i][r]
u[i*n+r]=a[i*n+r];
u[i*n+r]-=u[i*n+k]*u[k*n+r];
u[i*n+r]/=u[r*n+r];
//计算y，保存在x中
x[0]=b[0];
x[i]=b[i];
i;
x[i]-=u[i*n+k]*x[k];
//计算x，保存在x
x[n-1]/=u[（n-1）*n+n-1];
for（i=n-2;
i>
=0;
i--）
for（k=i+1;
x[i]/=u[i*n+i];
inti,j,n=3;
staticdoubleA[9]={1,2,3,2,5,2,3,1,5};
staticdoubleU[9];
staticdoubleX[3];
doubleB[3]={14,18,20};
triange（A,U,X,B,n）;
for（i=0;
for（j=0;
j<
j++）
printf（"
%g"
U[i*n+j]）;
x（%d）=%g\t"
i,X[i]）;
（2）高斯消去法
stdlib.h>
math.h>
intrgauss（intn,doublea[],doubleb[]）
boolflag=true;
//置非奇异标志
intk,i,j,is,p,q;
doubled,t;
int*js=（int*）malloc（n*sizeof（int））;
//记忆列交换信息
for（k=0;
k<
=n-2;
d=0.0;
for（i=k;
i<
=n-1;
i++）//全选主元
for（j=k;
j<
j++）
t=fabs（a[i*n+j]）;
if（t>
d）//记忆行、列交换信息
{
d=t;
js[k]=j;
is=i;
}
if（d+1.0==1.0）flag=0;
//置奇异标志
else
if（js[k]!
=k）
for（i=0;
i++）//列交换
p=i*n+k;
q=i*n+js[k];
t=a[p];
a[p]=a[q];
a[q]=t;
if（is!
for（j=k;
j<
j++）//行交换
p=k*n+j;
q=is*n+j;
t=b[k];
b[k]=b[is];
b[is]=t;
if（flag==0）//奇异返回
free（js）;
fail\n"
return0;
d=a[k*n+k];
for（j=k+1;
j++）//归一化
p=k*n+j;
a[p]=a[p]/d;
b[k]=b[k]/d;
for（i=k+1;
i++）//消元
for（j=k+1;
p=i*n+j;
a[p]=a[p]-a[i*n+k]*a[k*n+j];
b[i]=b[i]-a[i*n+k]*b[k];
d=a[（n-1）*n+n-1];
if（fabs（d）+1.0==1.0）//奇异返回
free（js）;
return0;
b[n-1]=b[n-1]/d;
//计算解向量的最后一个分量
i>
=0;
i--）//回代
t=0.0;
for（j=i+1;
t=t+a[i*n+j]*b[j];
b[i]=b[i]-t;
js[n-1]=n-1;
for（k=n-1;
k>
k--）//恢复解向量
if（js[k]!
=k）//解向量行交换
t=b[k];
b[k]=b[js[k]];
b[js[k]]=t;
free（js）;
return1;
//返回正常标志
//主函数程序
inti;
staticdoublea[9]={1,2,3,2,5,2,3,1,5};
staticdoubleb[3]={14,18,20};
if（rgauss（3,a,b）!
=0）
for（i=0;
3;
i++）
x（%d）=%g"
i,b[i]）;
3、实验结果
解的上面的矩阵存有L和U的信息。
从中可知L=
U=
（2）高斯消去法
实验总结：
1、这几次的实验使都是实现已有的算法，虽然不难但要用程序实现出来还是要对算法理解非常清楚，而且要很细心和足够的耐心去检查细小的错误。
2、通过实验常用的算法有了很深的认识，以后要主动去运用算法和基本思路解决问题。
升级会员 