2005年9月17日
#
可以仔细看看高斯消元法
方法1：
#include "stdio.h"
#include "stdlib.h"
void RKT(t,y,n,h,k,z)
int n; /*微分方程组中方程的个数，也是未知函数的个数*/
int k; /*积分的步数(包括起始点这一步)*/
double t; /*积分的起始点t0*/
double h; /*积分的步长*/
double y[]; /*存放n个未知函数在起始点t处的函数值,返回时,其初值在二维数组z的第零列中*/
double z[]; /*二维数组,体积为n x k.返回k个积分点上的n个未知函数值*/
{
extern void Func(); /*声明要求解的微分方程组*/
int i,j,l;
double a[4],*b,*d;
b=malloc(n*sizeof(double)); /*分配存储空间*/
if(b == NULL)
{
printf("内存分配失败\n");
exit(1);
}
d=malloc(n*sizeof(double)); /*分配存储空间*/
if(d == NULL)
{
printf("内存分配失败\n");
exit(1);
}
/*后面应用RK4公式中用到的系数*/
a[0]=h/2.0;
a[1]=h/2.0;
a[2]=h;
a[3]=h;
for(i=0; i<=n-1; i++)
z[i*k]=y[i]; /*将初值赋给数组z的相应位置*/
for(l=1; l<=k-1; l++)
{
Func(y,d);
for (i=0; i<=n-1; i++)
b[i]=y[i];
for (j=0; j<=2; j++)
{
for (i=0; i<=n-1; i++)
{
y[i]=z[i*k+l-1]+a[j]*d[i];
b[i]=b[i]+a[j+1]*d[i]/3.0;
}
Func(y,d);
}
for(i=0; i<=n-1; i++)
y[i]=b[i]+h*d[i]/6.0;
for(i=0; i<=n-1; i++)
z[i*k+l]=y[i];
t=t+h;
}
free(b); /*释放存储空间*/
free(d); /*释放存储空间*/
return;
}
main()
{
int i,j;
double t,h,y[3],z[3][11];
y[0]=-1.0;
y[1]=0.0;
y[2]=1.0;
t=0.0;
h=0.01;
RKT(t,y,3,h,11,z);
printf("\n");
for (i=0; i<=10; i++) /*打印输出结果*/
{
t=i*h;
printf("t=%5.2f\t ",t);
for (j=0; j<=2; j++)
printf("y(%d)=%e ",j,z[j][i]);
printf("\n");
}
}
void Func(y,d)
double y[],d[];
{
d[0]=y[1]; /*y0'=y1*/
d[1]=-y[0]; /*y1'=y0*/
d[2]=-y[2]; /*y2'=y2*/
return;
}
高斯消元法：
#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#include <math.h>
int GS(int,double**,double *,double);
double **TwoArrayAlloc(int,int);
void TwoArrayFree(double **);
void main()
{
int i,n;
double ep,**a,*b;
n = 3;
ep = 1e-4;
a = TwoArrayAlloc(n,n);
b = (double *)calloc(n,sizeof(double));
if(b == NULL)
{
printf("内存分配失败\n");
exit(1);
}
a[0][0]= 2; a[0][1]= 6; a[0][2]=-1;
a[1][0]= 5; a[1][1]=-1; a[1][2]= 2;
a[2][0]=-3; a[2][1]=-4; a[2][2]= 1;
b[0] = -12; b[1] = 29; b[2] = 5;
if(!GS(n,a,b,ep))
{
printf("不可以用高斯消去法求解\n");
exit(0);
}
printf("该方程组的解为：\n");
for(i=0;i<3;i++)
printf("x%d = %.2f\n",i,b[i]);
TwoArrayFree(a);
free(b);
}
int GS(n,a,b,ep)
int n;
double **a;
double *b;
double ep;
{
int i,j,k,l;
double t;
for(k=1;k<=n;k++)
{
for(l=k;l<=n;l++)
if(fabs(a[l-1][k-1])>ep)
break;
else if(l==n)
return(0);
if(l!=k)
{
for(j=k;j<=n;j++)
{
t = a[k-1][j-1];
a[k-1][j-1]=a[l-1][j-1];
a[l-1][j-1]=t;
}
t=b[k-1];
b[k-1]=b[l-1];
b[l-1]=t;
}
t=1/a[k-1][k-1];
for(j=k+1;j<=n;j++)
a[k-1][j-1]=t*a[k-1][j-1];
b[k-1]*=t;
for(i=k+1;i<=n;i++)
{
for(j=k+1;j<=n;j++)
a[i-1][j-1]-=a[i-1][k-1]*a[k-1][j-1];
b[i-1]-=a[i-1][k-1]*b[k-1];
}
}
for(i=n-1;i>=1;i--)
for(j=i+1;j<=n;j++)
b[i-1]-=a[i-1][j-1]*b[j-1];
return(1);
}
double **TwoArrayAlloc(int r,int c)
{
double *x,**y;
int n;
x=(double *)calloc(r*c,sizeof(double));
y=(double **)calloc(r,sizeof(double*));
for(n=0;n<=r-1;++n)
y[n]=&x[c*n];
return (y);
}
void TwoArrayFree(double **x)
{
free(x[0]);
free(x);
}
两个简单的曲线拟合的方法：
拉格朗日插值法：
#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
double LAG(int,double *,double *,double);
void main()
{
int n;
double *x,*y,t,lag;
t = 0.15;
n = 6;
x = (double*)calloc(n,sizeof(double));
if(x == NULL)
{
printf("内存分配失败\n");
exit(1);
}
y = (double*)calloc(n,sizeof(double));
if(y == NULL)
{
printf("内存分配失败\n");
exit(1);
}
x[0] = 0;
x[1] = 0.1;
x[2] = 0.195;
x[3] = 0.3;
x[4] = 0.401;
x[5] = 0.5;
y[0] = 0.39894;
y[1] = 0.39695;
y[2] = 0.39142;
y[3] = 0.38138;
y[4] = 0.36812;
y[5] = 0.35206;
lag = LAG(n,x,y,t);
printf("拉各朗日插值后得到的结果是：\n");
printf("f(%.2f)=%e\n",t,lag);
free(x);
free(y);
}
double LAG(n,x,y,t)
int n;
double *x;
double *y;
double t;
{
int i,j;
double p,s;
s = 0;
for(i=0;i<n-1;i++)
{
p = 1;
for(j=0;j<n-1;j++)
if(i!=j)
p*=(t-x[j])/(x[i]-x[j]);
s+=p*y[i];
}
return (s);
}
曲线拟合：
#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#include <math.h>
Smooth(double *,double *,double *,int,int,
double *,double *,double *);
void main()
{
int i,n,m;
double *x,*y,*a,dt1,dt2,dt3,b;
n = 20;
m = 6;
b = 0;
/*分别为x,y,a分配存贮空间*/
x = (double *)calloc(n,sizeof(double));
if(x == NULL)
{
printf("内存分配失败\n");
exit (0);
}
y = (double *)calloc(n,sizeof(double));
if(y == NULL)
{
printf("内存分配失败\n");
exit (0);
}
a = (double *)calloc(n,sizeof(double));
if(a == NULL)
{
printf("内存分配失败\n");
exit (0);
}
for(i=1;i<=n;i++)
{
x[i-1]=b+(i-1)*0.1;
/*每隔0.1取一个点，这样连续取n个点*/
y[i-1]=x[i-1]-exp(-x[i-1]);
/*计算x[i-1]点对应的y值作为拟合已知值*/
}
Smooth(x,y,a,n,m,&dt1,&dt2,&dt3); /*调用拟合函数*/
for(i=1;i<=m;i++)
printf("a[%d] = %.10f\n",(i-1),a[i-1]);
printf("拟合多项式与数据点偏差的平方和为：\n");
printf("%.10e\n",dt1);
printf("拟合多项式与数据点偏差的绝对值之和为：\n");
printf("%.10e\n",dt2);
printf("拟合多项式与数据点偏差的绝对值最大值为：\n");
printf("%.10e\n",dt3);
free(x); /*释放存储空间*/
free(y); /*释放存储空间*/
free(a); /*释放存储空间*/
}
Smooth(x,y,a,n,m,dt1,dt2,dt3 )
double *x; /*实型一维数组，输入参数，存放节点的xi值*/
double *y; /*实型一维数组，输入参数，存放节点的yi值*/
double *a; /*双精度实型一维数组，长度为m。返回m一1次拟合多项式的m个系数*/
int n; /*整型变量，输入参数，给定数据点的个数*/
int m; /*整型变量，输入参数，拟合多项式的项数*/
double *dt1; /*实型变量，输出参数，拟合多项式与数据点偏差的平方和*/
double *dt2; /*实型变量，输出参数，拟合多项式与数据点偏差的绝对值之和*/
double *dt3; /*实型变量，输出参数，拟合多项式与数据点偏差的绝对值最大值*/
{
int i,j,k;
double *s,*t,*b,z,d1,p,c,d2,g,q,dt;
/*分别为s,t,b分配存贮空间*/
s = (double *)calloc(n,sizeof(double));
if(s == NULL)
{
printf("内存分配失败\n");
exit (0);
}
t = (double *)calloc(n,sizeof(double));
if(t == NULL)
{
printf("内存分配失败\n");
exit (0);
}
b = (double *)calloc(n,sizeof(double));
if(b == NULL)
{
printf("内存分配失败\n");
exit (0);
}
z = 0;
for(i=1;i<=n;i++)
z=z+x[i-1]/n; /*z为各个x的平均值*/
b[0]=1;
d1=n;
p=0;
c=0;
for(i=1;i<=n;i++)
{
p=p+x[i-1]-z;
c=c+y[i-1];
}
c=c/d1;
p=p/d1;
a[0]=c*b[0];
if(m>1)
{
t[1]=1;
t[0]=-p;
d2=0;
c=0;
g=0;
for(i=1;i<=n;i++)
{
q=x[i-1]-z-p;
d2=d2+q*q;
c=y[i-1]*q+c;
g=(x[i-1]-z)*q*q+g;
}
c=c/d2;
p=g/d2;
q=d2/d1;
d1=d2;
a[1]=c*t[1];
a[0]=c*t[0]+a[0];
}
for(j=3;j<=m;j++)
{
s[j-1]=t[j-2];
s[j-2]=-p*t[j-2]+t[j-3];
if(j>=4)
for(k=j-2;k>=2;k--)
s[k-1]=-p*t[k-1]+t[k-2]-q*b[k-1];
s[0]=-p*t[0]-q*b[0];
d2=0;
c=0;
g=0;
for(i=1;i<=n;i++)
{
q=s[j-1];
for(k=j-1;k>=1;k--)
q=q*(x[i-1]-z)+s[k-1];
d2=d2+q*q;
c=y[i-1]*q+c;
g=(x[i-1]-z)*q*q+g;
}
c=c/d2;
p=g/d2;
q=d2/d1;
d1=d2;
a[j-1]=c*s[j-1];
t[j-1]=s[j-1];
for(k=j-1;k>=1;k--)
{
a[k-1]=c*s[k-1]+a[k-1];
b[k-1]=t[k-1];
t[k-1]=s[k-1];
}
}
*dt1=0;
*dt2=0;
*dt3=0;
for(i=1;i<=n;i++)
{
q=a[m-1];
for(k=m-1;k>=1;k--)
q=q*(x[i-1]-z)+a[k-1];
dt=q-y[i-1];
if(fabs(dt)>*dt3)
*dt3=fabs(dt);
*dt1=*dt1+dt*dt;
*dt2=*dt2+fabs(dt);
}
/*释放存储空间*/
free(s);
free(t);
free(b);
return(1);
}
如何求定积分
辛普生求积分：
#include <stdio.h>
#include <math.h>
double Function(double);
double SIMP1(double,double,int);
double SIMP2(double,double,double);
void main()
{
double a1,b1,eps;
int n1;
double Result1;
double Result2;
a1 = 0.0;
b1 = 0.8;
n1 = 4;
eps = 5e-7;
Result1 = SIMP1(a1,b1,n1);
Result2 = SIMP2(a1,b1,eps);
printf("利用定步长辛普生积分结果为：\n");
printf("I1 = %.10f\n",Result1);
printf("利用变步长辛普生积分结果为：\n");
printf("I2 = %e\n",Result2);
}
double SIMP1(a,b,n)
double a;
double b;
int n;
{
int i;
double h,s;
h=(a-b)/(2*n);
s=0.5*(Function(a)-Function(b));
for(i=1;i<=n;i++)
s+=2*Function(a+(2*i-1)*h)+Function(a+2*i*h);
return((b-a)*s/(3*n));
}
double SIMP2(a,b,eps)
double a;
double b;
double eps;
{
int k,n;
double h,t1,t2,s1,s2,p,x;
n=1;
h=b-a;
t1=h*(Function(a)+Function(b))/2;
s1 = t1;
while(1)
{
p = 0;
for(k=0;k<=n;k++)
{
x = a+(k+0.5)*h;
p+=Function(x);
}
t2=(t1+h*p)/2;
s2=(4*t2-t1)/3;
if(fabs(s2-s1)>=eps)
{
t1=t2;
n=n+n;
h=h/2;
s1=s2;
continue;
}
break;
}
return(s2);
}
double Function(double x)
{
return(cos(x));
}
欧拉方程法：
#include "stdio.h"
#include "stdlib.h"
#include <math.h>
int Func(y,d)
double y[],d[];
{
d[0]=y[1]; /*y0'=y1*/
d[1]=-y[0]; /*y1'=y0*/
d[2]=-y[2]; /*y2'=y2*/
return(1);
}
void Euler1(t,y,n,h,k,z)
int n; /*整型变量，微分方程组中方程的个数，也是未知函数的个数*/
int k; /*整型变量。积分步数(包括起始点这一步)*/
double t; /*双精度实型变量,对微分方程进行积分的起始点t0*/
double h; /*双精度实型变量。积分步长*/
double y[]; /*双精度实型一维数组，长度为n。存放n个未知函数yi在起始点t0处的函数值*/
double z[]; /*双精度实型二维数组，体积为nxk。返回k个积分点(包括起始点)上的未知函数值*/
{
extern int Func();
int i,j;
double *d;
d=malloc(n*sizeof(double));
if(d == NULL)
{
printf("内存分配失败\n");
exit(1);
}
/*将方程组的初值赋给数组z[i*k]*/
for (i=0; i<=n-1; i++)
z[i*k]=y[i];
for (j=1; j<=k-1; j++)
{
Func(y,d); /*求出f(x)*/
for(i=0; i<=n-1; i++)
y[i]=z[i*k+j-1]+h*d[i];