yingkou的Blog

求解微分方程

yingkou的Blog — Sat, 17 Sep 2005 05:23:00 GMT

可以仔细看看高斯消元法

方法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);

}

曲线拟合

yingkou的Blog — Sat, 17 Sep 2005 05:21:00 GMT

两个简单的曲线拟合的方法：

拉格朗日插值法：

#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);

}