全区间积分的双边法(常微分方程组的求解)
来源:互联网 发布:蓝桥杯java b组怎么样 编辑:程序博客网 时间:2024/05/17 05:34
/*
代码作者:不详
代码整理者:设计天下 MySDN网站 算法天下工作室
功能:全区间积分的双边法(常微分方程组的求解)
*/
#include "stdio.h"
#include "stdlib.h"
#include "math.h"
/*全区间积分的定步长欧拉方法*/
/*全区间积分的维梯方法*/
/*全区间积分的定步长龙格-库塔方法*/
typedef struct _fode {
int n; /*微分方程组中方程个数,也是未知函数的个数*/
int steps; /*积分步数(包括起始点这一步)*/
double lens; /*积分的步长*/
double t; /*对微分方程进行积分的起始点 */
double *y; /*存放n个未知函数在起始点t处的函数值*/
double *z; /*返回steps个积分点(包括起始点)上的未知函数值*/
void (*ptr)(); /*指向计算微分方程组中各方程右端函数值的函数名(由用户自编)*/
} FODE, *FODEP;
/*全区间积分的变步长默森方法*/
/*全区间积分的双边法*/
/*全区间积分的阿当姆斯预报校正法*/
/*全区间积分的哈明方法*/
typedef struct _fode2 {
int n; /*微分方程组中方程个数,也是未知函数的个数*/
int steps; /*积分步数(包括起始点这一步)*/
double lens; /*积分的步长*/
double eps; /*控制精度要求*/
double t; /*对微分方程进行积分的起始点 */
double *y; /*存放n个未知函数在起始点t处的函数值*/
double *z; /*返回steps个积分点(包括起始点)上的未知函数值*/
void (*ptr)(); /*指向计算微分方程组中各方程右端函数值的函数名(由用户自编)*/
} FODE2, *FODE2P;
/*积分一步的特雷纳方法*/
typedef struct _tlode {
int n; /*微分方程组中方程个数,也是未知函数的个数*/
double lens; /*积分的步长*/
double t; /*对微分方程进行积分的起始点 */
double *y; /*存放n个未知函数在起始点t处的函数值。返回t+lens点处的n个未知函数值*/
void (*ptr)(); /*指向计算微分方程组中各方程右端函数值的函数名(由用户自编)*/
} TODE, *TODEP;
/*积分一步的变步长欧拉方法*/
/*积分一步的变步长龙格-库塔方法*/
/*积分一步的连分式法*/
typedef struct _eode {
int n; /*微分方程组中方程个数,也是未知函数的个数*/
double lens; /*积分的步长*/
double eps; /*积分的精度要求*/
double t; /*对微分方程进行积分的起始点 */
double *y; /*存放n个未知函数在起始点t处的函数值。返回t+lens点处的n个未知函数值*/
void (*ptr)(); /*指向计算微分方程组中各方程右端函数值的函数名(由用户自编)*/
} EODE, *EODEP;
/*积分一步的变步长基尔方法*/
typedef struct _eode2 {
int n; /*微分方程组中方程个数,也是未知函数的个数*/
double lens; /*积分的步长*/
double eps; /*积分的精度要求*/
double t; /*对微分方程进行积分的起始点 */
double *y; /*存放n个未知函数在起始点t处的函数值。返回t+lens点处的n个未知函数值*/
double *q; /*在主函数第一次调用本函数时,应赋值0,
以后每调用一次函数(即每积分一步),将由本函数的返回值以便循环使用。*/
void (*ptr)(); /*指向计算微分方程组中各方程右端函数值的函数名(由用户自编)*/
} EODE2, *EODE2P;
/*二阶微分方程边值问题的数值解法*/
typedef struct _dode {
int n; /*求解区间[a,b]的等分点数(包括左端点a与右端点b)*/
double a; /*求解区间的左端点*/
double b; /*求解区间的右端点求*/
double ya; /*未知函数在求解区间左端点处的函数值y(a)*/
double yb; /*未知函数在求解区间右端点处的函数值y(b)*/
double *y; /*返回n个等距离散点上的未知函数值*/
void (*ptr)(); /*指向计算二阶微分方程中函数值的函数名(由用户自编)*/
} DODE, *DODEP;
static void bilateral_method_lis(int n,double eps,double t,double lens,double y[],void (*ptr)())
{
int m,i,j,k;
double hh,p,dt,x,tt,q,a[4];
double *g,*b,*c,*d,*e;
g=malloc(n*sizeof(double));
b=malloc(n*sizeof(double));
c=malloc(n*sizeof(double));
d=malloc(n*sizeof(double));
e=malloc(n*sizeof(double));
hh=lens;
m=1;
x=t;
for (i=0; i<n; i++)
{
c[i]=y[i];
}
do
{
a[1]=a[0]=hh/2.0;
a[3]=a[2]=hh;
for (i=0; i<n; i++)
{
g[i]=y[i];
y[i]=c[i];
}
dt=lens/m; t=x;
for (j=0; j<m; j++)
{
(*ptr)(t,y,n,d);
for (i=0; i<n; i++)
{
b[i]=y[i]; e[i]=y[i];
}
for (k=0; k<3; k++)
{
for (i=0; i<n; i++)
{
y[i]=e[i]+a[k]*d[i];
b[i]+=a[k+1]*d[i]/3.0;
}
tt=t+a[k];
(*ptr)(tt,y,n,d);
}
for (i=0; i<n; i++)
{
y[i]=b[i]+hh*d[i]/6.0;
}
t+=dt;
}
p=0.0;
for (i=0; i<n; i++)
{
q=fabs(y[i]-g[i]);
if (q>p)
{
p=q;
}
}
hh=hh/2.0;
m<<=1;
}
while (p>=eps);
/*释放动态分配的内存*/
free(g); free(b); free(c); free(d); free(e);
return;
}
void bilateral_lis(FODE2P ap)
{
int n,steps,i,j;
double lens,eps,t,t0,qq;
double *d,*p,*u,*v,*w,*y,*z;
n=ap->n;
steps=ap->steps;
d=malloc(n*sizeof(double));
p=malloc(n*sizeof(double));
u=malloc(n*sizeof(double));
v=malloc(n*sizeof(double));
w=malloc(n*sizeof(double));
lens=ap->lens;
eps=ap->eps;
y=ap->y;
z=ap->z;
for (i=0; i<n; i++)
{
p[i]=0.0;
z[i*steps]=y[i];
}
t0=t=ap->t;
(*ap->ptr)(t,y,n,d);
for (j=0; j<n; j++)
{
u[j]=d[j];
}
bilateral_method_lis(n,eps,t,lens,y,ap->ptr);
t+=lens;
(*ap->ptr)(t,y,n,d);
for (j=0; j<n; j++)
{
z[j*steps+1]=y[j];
v[j]=d[j];
}
for (j=0; j<n; j++)
{
p[j]=-4.0*z[j*steps+1]+5.0*z[j*steps]+2.0*lens*(2.0*v[j]+u[j]);
y[j]=p[j];
}
t+=lens;
(*ap->ptr)(t,y,n,d);
for (j=0; j<n; j++)
{
qq=2.0*lens*(d[j]-2.0*v[j]-2.0*u[j])/3.0
+4.0*z[j*steps+1]-3.0*z[j*steps];
z[j*steps+2]=(p[j]+qq)/2.0;
y[j]=z[j*steps+2];
}
for (i=3; i<steps; i++)
{
t=t0+(i-1)*lens;
(*ap->ptr)(t,y,n,d);
for (j=0; j<=n-1; j++)
{
u[j]=v[j];
v[j]=d[j];
}
for (j=0; j<n; j++)
{
qq=-4.0*z[j*steps+i-1]+5.0*z[j*steps+i-2];
p[j]=qq+2.0*lens*(2.0*v[j]+u[j]);
y[j]=p[j];
}
t+=lens;
(*ap->ptr)(t,y,n,d);
for (j=0; j<n; j++)
{
qq=2.0*lens*(d[j]-2.0*v[j]-2.0*u[j])/3.0
+4.0*z[j*steps+i-1]-3.0*z[j*steps+i-2];
y[j]=(p[j]+qq)/2.0;
z[j*steps+i]=y[j];
}
}
/*释放动态分配的内存*/
free(d);
free(p);
free(u);
free(v);
free(w);
return;
}
void bilateral_ptr_lis(double t,double y[],int n,double d[])
{
t=t;
n=n;
d[0]=-2.0*y[1];
d[1]=2.0*y[0];
return;
}
int main()
{
double y[2]={0.0, 1.0};
double z[2][11]={0};
FODE2 fa={2, 11, 0.1, 0.00001, 0.1, y, (double*)z, bilateral_ptr_lis};
int i;
bilateral_lis(&fa);
printf("/n");
for (i=0; i<fa.steps; i++)
{
printf("t=%7.3f y(0)=%e y(1)=%e/n", i*fa.lens,z[0][i],z[1][i]);
}
printf("/n");
return 0;
}
- 全区间积分的双边法(常微分方程组的求解)
- 全区间积分的阿当姆斯预报校正法(常微分方程组的求解)
- 全区间积分的哈明方法(常微分方程组的求解)
- 全区间积分的定步长欧拉方法(常微分方程组的求解)
- 求解一个有趣的常微分方程组
- 求解微分方程组的ODE算法
- Matlab求解微分方程组
- 解常微分方程组 python
- Matlab求解中性类型的时滞微分方程组-中性类型的时滞微分方程
- 用python实现解常微分方程组的简单示例以及用odeint解常微分方程的范例
- matlab求解常微分方程组/传染病模型并绘制SIR曲线
- 线性微分方程组的公式解法
- Rung-Kutta解线性常微分方程组
- 用显式欧拉格式和改进的欧拉格式求解常微分初值问题
- 多元一次方程组的求解
- ADC的积分非线性和微分非线性
- ADC的积分非线性和微分非线性
- 微分 积分 时间常数的 限制 和图形
- 多元线性回归分析
- 全区间积分的哈明方法(常微分方程组的求解)
- asp.net 创建和终止线程(多线程)
- 全区间积分的阿当姆斯预报校正法(常微分方程组的求解)
- 当前全量表与拉链表关联查询获取历史时刻准确数据实例
- 全区间积分的双边法(常微分方程组的求解)
- Windows系统密码破解全攻略(hash破解)
- 对销毁对话框的分析
- 全区间积分的定步长欧拉方法(常微分方程组的求解)
- 简明x86汇编语言教程(1)
- 今日任务-调度航信预警规则+调度操作单的设计
- 利用MIXER获取麦克风录音音量
- 随机样本分析
- QListWidget的edit如何及时响应?(thinkvd开发日志)