解微分方程算法C语言程序（精心整理版）
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<graphics.h>
float k1,k2,k3,k4,m1,m2,m3,m4,n1,n2,n3,n4;
float t,ts,tf,h;
float x[1001],y[1001],z[1001],x_axis[1001];
int count;
void initial( )
{
x[0] = 1.0 ;
y[0] = -1.0;
z[0] = 0.0;
ts = 0.0;
tf = 20.0;
count = 500;
}
void runge( )
{
int i = 0;
h = (tf - ts)/count;
for(i = 1;i <= count;i++)
{
t = ts + (i-1)*h;
k1 = -x[i-1] + 2*y[i-1] + 6*z[i-1];
m1 = -y[i-1] + 3*z[i-1] + 2*sin(t);
n1 = -z[i-1] + sqrt(t)*exp(-t) + cos(t);
k2 = -(x[i-1] + k1*(h/2.0)) + 2*(y[i-1] + m1*(h/2.0)) + 6*(z[i-1] + n1*(h/2.0));
m2 = -(y[i-1] + m1*(h/2.0)) + 3*(z[i-1] + n1*(h/2.0)) + 2*sin(t + h/2.0);
n2 = -(z[i-1] + n1*(h/2.0)) + sqrt(t + h/2.0)*exp(-(t + h/2.0)) + cos(t + h/2.0);
k3 = -(x[i-1] + k2*(h/2.0)) + 2*(y[i-1] + m2*(h/2.0)) + 6*(z[i-1] + n2*(h/2.0));
m3 = -(y[i-1] + m2*(h/2.0)) + 3*(z[i-1] + n2*(h/2.0)) + 2*sin(t + h/2.0);
n3 = -(z[i-1] + n2*(h/2.0)) + sqrt(t + h/2.0)*exp(-(t + h/2.0)) + cos(t + h/2.0);
k4 = -(x[i-1] + k3*(h/2.0)) + 2*(y[i-1] + m3*(h/2.0)) + 6*(z[i-1] + n3*(h/2.0));
m4 = -(y[i-1] + m3*(h/2.0)) + 3*(z[i-1] + n3*(h/2.0)) + 2*sin(t + h/2.0);
n4 = -(z[i-1] + n3*(h/2.0)) + sqrt(t + h/2.0)*exp(-(t + h/2.0)) + cos(t + h/2.0);
x[i] = x[i-1] + h*(k1 + 2*k2 + 2*k3 + k4)/6.0;
y[i] = y[i-1] + h*(m1 + 2*m2 + 2*m3 + m4)/6.0;
z[i] = z[i-1] + h*(n1 + 2*n2 + 2*n3 + n4)/6.0;
}
明解c语言
}
int main( )
{
int i,j,k,gdriver,gmode;
char s[30];
initial( );
runge( );
for(i = 0;i <= count;i++)
x_axis[i] = i;
gdriver=DETECT; /*图形显示初始化*/
initgraph(&gdriver,&gmode,"c:\\tc");
setbkcolor(9);
outtextxy(240, 40, "Runge Rutta Arithmetic"); /*标题*/
/*outtextxy(400, 445, "THE PRODUCER WANGHUI "); 制作人*/
outtextxy(150, 415, "white is x red is y blue is z");
rectangle(20, 70, 620, 430); /*各种矩形框*/
rectangle(20, 70, 620, 410);
line(25, 250, 595, 250); /*建立X1轴坐标*/
line(595, 250, 590, 248);
line(595, 250, 590, 252);
outtextxy(580, 255, "t/s");
/*line(330, 405, 330, 75); 建立Y1轴坐标
line(330, 75, 328, 80);
line(330, 75, 332, 80);
outtextxy(335, 80, "h/mm");*/
line(40,405,40,75); /*建立Y1轴坐标*/
line(40,75,38,80);
line(40,75,42,80);
outtextxy(45, 80, "h/mm");
/*sprintf(s, "%d", x[0]);
outtextxy(18,445, s);
sprintf(s, "%d", x[1]);
outtextxy(26,445, s);
sprintf(s, "%d", x[2]);
outtextxy(34,445, s);
sprintf(s, "%d", x[3]);
outtextxy(42,445, s);*/
for(i = 0;i <= count;i++)
{
putpixel(40 + x_axis[i],250 - 15*x[i],15);
putpixel(40 + x_axis[i],250 - 15*y[i],4);
putpixel(40 + x_axis[i],250 - 15*z[i],1);
setcolor(15);
delay(1000); /*延时10000ms,动态输出*/
}
getch();
return(0);
}

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。