Question

求 解三阶微分方程的数值方法

三阶常微分方程
0.001y'''+0.11y''+y'+10y=10

求数值解法，如欧拉，四阶龙格库塔等，精度越高越好。
因为需要用C编程，所以尽可能给出具体公式。

如果有程序，给加分
初始值 x=0 y=0

Answer 1

初始值，得给一个。
一般一元多阶微分方程要化成：多元一阶微分方程组。

/*===================================================
Author : Wacs5
Date : 2010-01-05(YYYY-MM-DD)
Function: 四阶龙格库塔法
Compiler: TC2.0
Output : 生成rk4_data.txt文件
==================================================*/
#include <math.h>
#include <stdio.h>
#include <conio.h>
#include <stdlib.h>
#define EQUN 3

main()
{
int i,j;
int n=1000;
double t0=0,tend=1;
double y0[EQUN]={0.99,0,0}; /*初值*/

int fun_dequ(double t,double *y,int n,double *k);
int RK4_vn(double t0,double tend,int t_n,int equ_n,double *y0);
int RK4_vn_step(double t,double dt,int n,double *y0,double *y1);

system("cls");
RK4_vn(t0,tend,n,EQUN,y0);
getch();
}

/*===================================================================
定步长四阶龙格求解微分方程
==================================================================*/
int RK4_vn(double t0,double tend,int t_n,int equ_n,double *y0)
{
int i,j;
double dt=(tend-t0)/t_n;
double t;
double *py1,*pyhead;
int fun_dequ(double t,double *y,int n,double *k);
int RK4_vn_step(double t,double dt,int n,double *y0,double *y1);
FILE *fp;

py1=(double *)malloc(equ_n*(t_n+1)*sizeof(double));
if (py1==NULL)
{
printf("Malloc error:\nPress any key to exit:");
getch();
exit(1);
}
pyhead=py1; /*backup*/

for (i=0;i<equ_n;i++)
py1[i]=y0[i];

for (i=0,t=t0;i<t_n;i++,py1+=equ_n,t+=dt)
RK4_vn_step(t,dt,equ_n,py1,py1+equ_n);

/*save the data*/
fp=fopen("RK4_data.txt","w");
if (fp==NULL)
{
printf("File Open Error:\nPress any key to exit:");
getch();
exit(1);
}
for (i=0,py1=pyhead,t=t0;i<=t_n;i++,t+=dt,fprintf(fp,"\n"))
for (j=0,fprintf(fp,"%g\t",t);j<equ_n;j++,fprintf(fp,"\t"))
fprintf(fp,"%g",*py1++);
fclose(fp);

py1-=3;
printf("t=%f\n",tend);
for (i=0;i<EQUN;i++)
printf("y[%d]=%f\n",i,*py1++);
printf("OK");

return 0;
}

/*===================================================================
定步长四阶龙格中的一个时间步长
==================================================================*/
int RK4_vn_step(double t,double dt,int equ_n,double *y0,double *y1)
{
int i,j;
static int first=0;
double a[5]={0,0.5,0.5,1,0};
double b[5]={0,1.0/6,1.0/3,1.0/3,1.0/6};
static double *ki,*ymid;
int fun_dequ(double t,double *y,int n,double *k);

if (first==0)
{
first=1;
ki=(double *)malloc(equ_n*sizeof(double));
ymid=(double *)malloc(equ_n*sizeof(double));
if (ki==NULL || ymid==NULL)
{
printf("Malloc error:\nPress any key to exit:");
getch();
exit(1);
}
}
for (i=0;i<equ_n;i++)
{
ki[i]=0;
y1[i]=y0[i];
}

for (i=0;i<5;i++)
{
a[i]*=dt;
b[i]*=dt;
}

for (i=0;i<4;i++)
{
for (j=0;j<equ_n;j++)
ymid[j]=y0[j]+a[i]*ki[j];
fun_dequ(t+a[i],ymid,equ_n,ki);
for (j=0;j<equ_n;j++)
y1[j]+=b[i+1]*ki[j];
}
return 1;
}

/*===================================================================
描述微分方程的右端项f(t,yn)
==================================================================*/
int fun_dequ(double t,double *y,int equ_n,double *k)
{
int i;
k[0]=y[1];
k[1]=y[2];
k[2]=10000-110*y[2]-1000*y[1]-10000*y[0];
return 1;
}

Answer 2

还好是三阶的...
先来看x^3+px+q=0
由卡尔丹公式：
x1=(-q/2+((q/2)^2+(p/3)^3)^(1/2))^(1/3)+(-q/2-((q/2)^2+(p/3)^3)^(1/2))^(1/3)
x2=w(-q/2+((q/2)^2+(p/3)^3)^(1/2))^(1/3)+w^2(-q/2-((q/2)^2+(p/3)^3)^(1/2))^(1/3)
x3=w^2(-q/2+((q/2)^2+(p/3)^3)^(1/2))^(1/3)+w(-q/2-((q/2)^2+(p/3)^3)^(1/2))^(1/3)
其中 w=(-1+i3^(1/2))/2,w^2=(-1-i3^(1/2))/2

再来看ax^3+bx^2+cx+d=0
上式除以a并设x=y-b/3a,转化成 y^3+py+1=0的形式，求出y1,y2,y3后有
x1=y1-b/3a,x2=y2-b/3a,x2=y1-b/3a

这样就得到了特征根，构造通解后再求特解应该没有障碍，不过编程的话，量可不小，考虑用下成熟的数学包吧。

Answer 3

由卡尔丹公式：
x1=(-q/2+((q/2)^2+(p/3)^3)^(1/2))^(1/3)+(-q/2-((q/2)^2+(p/3)^3)^(1/2))^(1/3)
x2=w(-q/2+((q/2)^2+(p/3)^3)^(1/2))^(1/3)+w^2(-q/2-((q/2)^2+(p/3)^3)^(1/2))^(1/3)
x3=w^2(-q/2+((q/2)^2+(p/3)^3)^(1/2))^(1/3)+w(-q/2-((q/2)^2+(p/3)^3)^(1/2))^(1/3)
其中 w=(-1+i3^(1/2))/2,w^2=(-1-i3^(1/2))/2

由ax^3+bx^2+cx+d=0可知
上式除以a并设x=y-b/3a,转化成 y^3+py+1=0的形式，求出y1,y2,y3后有
x1=y1-b/3a,x2=y2-b/3a,x2=y1-b/3a

即得特征根
名词解释：
微分方程指描述未知函数的导数与自变量之间的关系的方程。微分方程的解是一个符合方程的函数。而在初等数学的代数方程，其解是常数值。
微分方程的应用十分广泛，可以解决许多与导数有关的问题[1]:p.1。物理中许多涉及变力的运动学、动力学问题，如空气的阻力为速度函数的落体运动等问题，很多可以用微分方程求解。此外，微分方程在化学、工程学、经济学和人口统计等领域都有应用。
数学领域对微分方程的研究着重在几个不同的面向，但大多数都是关心微分方程的解。只有少数简单的微分方程可以求得解析解。不过即使没有找到其解析解，仍然可以确认其解的部份性质。在无法求得解析解时，可以利用数值分析的方式，利用电脑来找到其数值解。 动力系统理论强调对于微分方程系统的量化分析，而许多数值方法可以计算微分方程的数值解，且有一定的准确度。