2019数值分析实验报告存档

一、插值与拟合

#include<stdio.h>
#include<stdlib.h>
double Newton(double *,double *,double );
void main (void)
{
	double x[3]={1.000000,4.000000,9.000000};
	double y[3]={1.000000,2.000000,3.000000};
	double u;
	printf("请输入预测点x: ");
	scanf("%lf",&u);
	printf("预测点的近似函数值为：%lf",Newton(x,y,u)); 
}
double Newton(double *a,double *b,double p)
{
	double add=0,t1=1;
	double diQuo[2];//差商 
	int i,j,k;
	for(k=0;k<2;k++)
	{
		add=0;
	   	for(i=0;i<=k+1;i++)
      	{
      		t1=1;
	        j=0;
	    	while(j<=k+1)
	    	{
    			if(j!=i)
      			{
     				t1*=(a[i]-a[j]);
	     			j++;
	    		}
	    		else
    			{
	    			j++;
	    			continue;
	    		}
    		}
     		add+=(b[i]/t1);
     	    printf("%lf,",t1);
    	}
    	diQuo[k]=add;
	}
	add=b[0]+diQuo[0]*(p-a[0])+diQuo[1]*(p-a[0])*(p-a[1]);
	return add;
	
}

二、数值积分

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
double Trapezoid(double,double,int);//复化梯形公式 
double Simpson(double ,double ,int);//复化辛普生公式 
double Romberg(double );//龙贝格算法 
void main(void)
{
	double a,b,r;
	int n;
	printf("-----------------------------复化梯形公式-------------------------\n"); 
	printf("此程序使用的函数为f(x)=sinx/x"); 
	printf("请输入积分下限a:");
	scanf("%lf",&a);
	printf("\n请输入积分上限b:");
	scanf("%lf",&b);
	printf("请输入划分：");
	scanf("%d",&n);
	printf("复化梯形公式求得函数在%f到%f上的积分值为:%lf\n",a,b,Trapezoid(a,b,n));
	printf("-----------------------------复化辛普生公式-----------------------\n");
	printf("此程序使用的函数为f(x)=sinx/x"); 
	printf("请输入积分下限a:");
	scanf("%lf",&a);
	printf("\n请输入积分上限b:");
	scanf("%lf",&b);
	printf("请输入划分n:");
	scanf("%d",&n);
	printf("复化辛普生公式求得函数在%f到%f上的积分值为:%lf\n",a,b,Simpson(a,b,n));
	printf("-----------------------------龙贝格算法---------------------------\n");
	printf("此程序使用的函数为f(x)=sinx/x"); 
	printf("请输入误差限：");
	scanf("%lf",&r);
    printf("龙贝格算法求得积分值为%lf",Romberg(r));
 } 
double Trapezoid(double a1,double b1,int n2)
{
	double T,x,h,add,fa,fb;
	int i=0;
	add=0;
	h=(b1-a1)/n2;
	fa=1;
	fb=sin(b1)/b1;
	for(i=1;i<=n2-1;i++)
	{
		x=a1+i*h;
		add=add+(sin(x)/x);
	}
	printf("%lf",add);
	T=(h/2)*(fa+2*add+fb);
	return T;
}
double Simpson(double a2,double b2,int n1)
{
	double h,x,add,fa,fb,f1,f2,f3;
	int i;
	h=(b2-a2)/n1;
	add=(h/6)*(1+4*sin(a2+h/2)/(a2+h/2)+sin(a2+h)/(a2+h));
	x=a2+h;
	for(i=1;i<n1;i++)
	{
		f1=sin(x)/x;
		x=x+h/2;
		f2=sin(x)/x;
		x=x+h/2;
		f3=sin(x)/x;
		add+=(h/6)*(f1+4*f2+f3);
	}
	return add;
}
double Romberg(double r1)
{
	double x,h,S,S1,S2,T1,T2,C1,C2,R1,R2;
	int k;
	h=1;
	T1=1;
	k=1;
	do
	{
		S=0;
		x=0+h/2;
		do
		{
			S=S+(sin(x)/x);
			x=x+h;
		}while(x<1);
		T2=T1/2+h*S/2;
		S2=T2+(T2-T1)/3;
		if(k==1)
		{
			k++;
			h=h/2;
			T1=T2;
			S1=S2;
			S=0;
		    x=0+h/2;
		    do
		    {
			    S=S+(sin(x)/x);
			    x=x+h;
		    }while(x<1);
		    T2=T1/2+h*S/2;
		    S2=T2+(T2-T1)/3;
		}
		C2=S2+(S2-S1)/15;
		if(k==2)
		{
			C1=C2;
			k++;
			h=h/2;
			T1=T2;
			S1=S2;
			S=0;
		    x=0+h/2;
		    do
		    {
			    S=S+(sin(x)/x);
			    x=x+h;
		    }while(x<1);
		    T2=T1/2+h*S/2;
		    S2=T2+(T2-T1)/3;
		    C2=S2+(S2-S1)/15;
		}
		R2=C2+(C2-C1)/63;
	    if(k==3)
		{
			R1=R2;
			C1=C2;
			k++;
			h=h/2;
			T1=T2;
			S1=S2;
			S=0;
		    x=0+h/2;
		    do
		    {
			    S=S+(sin(x)/x);
			    x=x+h;
		    }while(x<1);
		    T2=T1/2+h*S/2;
		    S2=T2+(T2-T1)/3;
		    C2=S2+(S2-S1)/15;
		    R2=C2+(C2-C1)/63;
		} 
	}while(fabs(R2-R1)>=r1);
	return R2;
}

三、非线性方程求根迭代法

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
double fun1(double x)//定义原函数 
{
    double y;
	y=pow(x,3)-x-1;
	return y;	
} 
double fun2(double x)//定义导函数 
{
	double y;
	y=3*x*x-1;
	return y;
}
void Newton(double x,double R,double N,double M)
{
	int k=0;
	int i;
	double n;
	double x1;
	while(k<N)
	{
		if(fun2(x)==0)
    	{
	    	printf("奇异标志\n");
    		return;
	    }
    	else
    	{
     		i=0;
    		n=1;
    		while(1)
    		{
    			x1=x-n*fun1(x)/fun2(x);
     			printf("%10f\n",x1);
    			if(fabs( fun1(x1) ) < fabs( fun1(x) ) )
     			{
    				break;
    			}
    			else
    			{
    				i++;
    				n=n/2;
     			}
     			if(i>=M)
      			{
    				printf("重新输入迭代初值\n");
     				return;
    			}
	    	}  
		    if(fabs(x1-x)<R)
    		{
    			printf("结果为:%10f\n",x1);
                return;
    		}
    		else
    		{
    			k++;
    			x=x1;
    		}
    	}
	}
	printf("迭代失败\n");
	return;
 } 
void main(main)
{
	double x,r;
	int N1,N2;
	printf("请输入迭代初值：");
	scanf("%f",&x);
	printf("请输入误差限:");
	scanf("%f",&r);
	printf("请输入迭代最大次数：");
	scanf("%d",&N1);
	printf("请输入下山最大次数：");
	scanf("%d",&N2);
	Newton(x,r,N1,N2); 
}
 
 
 
 
 
 
 
 

四、数值微分

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
double Euler(double ,double ,double );
double Euler2(double ,double ,double );
double RungKutta(double ,double ,double );
void main(void)
{
	int i;
    double h=0.1;
    double y1,y2,y3;
	y1=y2=y3=1;
	double x=0.000000;
	printf("y'=y-2x/y   0<x<1\n");
	printf("y(0)=1\n");
	printf("步长h=0.1\n");
	printf("   x        y(四阶龙格) y(欧拉) y(改进欧拉) y(精确)        \n");
	for(i=0;i<10;i++)
	{
		y1=RungKutta(h,x,y1);
		y2=Euler(h,x,y2);
		y3=Euler2(h,x,y3);
		printf("%10f",x+h);
		printf("%10f",y1);
		printf("%10f",y2);
		printf("%10f",y3);
		printf("%10f\n",sqrt(1+2*(x+h)));
		x+=0.100000;
	}
	
}
double Euler(double h,double x0,double y0)
{
	double yn;
	yn=y0+h*(y0-(2*x0/y0));
	return yn;
}
double Euler2(double h,double x0,double y0)
{
	double yp,yc,yn,xn;
	yp=y0+h*(y0-2*x0/y0);
	xn=x0+h;
	yc=y0+h*(yp-2*xn/yp);
	yn=(yp+yc)/2;
	return yn;
}
double RungKutta(double h,double x0,double y0)
{
	double K1=0,K2=0,K3=0,K4=0;
    double yn;
	K1=y0-2*x0/y0;
	K2=y0+h*K1/2-(2*x0+h)/(y0+h*K1/2);
	K3=y0+h*K2/2-(2*x0+h)/(y0+h*K2/2);
	K4=y0+h*K3-2*(x0+h)/(y0+h*K3);
	yn=y0+h*(K1+2*K2+2*K3+K4)/6;
	return yn;
}