龙哥库他发matlab程序,编程实现四阶龙哥库塔法解方程

% program's name :Solving nonhomogeneous equation

% designer       : Lv Shiqi

% date           :2014-03-14

%创建了一个10*10的矩阵，里面的数据是以1为公差的等差数列，从1一直加到99

y=0:99;

%这里的h和k都可以改变。图形就可以进行相应改变。

% k是输入函数的系数，h是时间间隔

h=0.1;

k=2;

dis=h/2.0;

pre=h/6.0;

%和C语言不同，matlab的数组下标是从1开始的。

for i=1:99

k1=k-y(i);

k2=k-(y(i)+dis*k1);

k3=k-(y(i)+dis*k1);

k4=k-(y(i)+h*k3);

y(i+1)=y(i)+(k1+2*k2+2*k3+k4)*pre;

end

%这里表示x是以0.2为公差的等差数列。从0一直加到19.8，一共100个数

x=0:0.1:9.9;

%绘图函数plot，只需要将x与y值传入其中即可绘图。

plot(x,y)

二、C源程序

// Program's name :Solving nonhomogeneous equation

// Designer       :Lv Shiqi

// Date           :2014-03-14

// Version    :Ver0.2

#include

#define FAILED    -1

// the time interval is 0.1s, it can be changed

#define h  0.1

// the parameters of the output functiong is 2, it canbe changed

#define k  2

int main()

{

int i;

double k1, k2, k3, k4;

double y[100] ={0.0};

double dis =h/2.0;

double pre =h/6.0;

FILE   *fp;

// ifdata.txt does not exist, exit the program

if((fp=fopen("data.txt","w"))==NULL){

printf("Cannotopen the file!\n");

returnFAILED;

}

for(i = 0;i<99; i++){

k1  = k-y[i];

k2  = k-(y[i]+dis*k1);

k3  = k-(y[i]+dis*k2);

k4  = k-(y[i]+h*k3);

y[i+1] = y[i]+(k1+2*k2+2*k3+k4)*pre;  // 对应Yn+1=Yn+(K1+2K2+2K3+K4)*(h/6)

}

for(i = 0;i<100; i++){

fprintf(fp,"%f\t%f\n",i*h, y[i]);

}

fclose(fp);

system("./mysql.sh");

return 0;

}

三、matlab实现的图形

四、C语言实现的图形

五、数据分析

本次程序设计中求解的方程式为：Ty’(t)+y(t)=Ku(t)，其中u(t)是阶跃函数。在我的程序中T=1，K=2，y(0)=0。

通过上部分展现的图形可以看出，图形最终会趋近于一个常数，而这个常数就是K，即2，而且，图像增长的速度是先快后慢。而且，在测试程序的时候我改小了T值，图像会更陡，改大了T值，图像会更加平缓。

而通过人工计算，方程解应该是y(t)=K(1-et/T)，通过分析，这个函数的图像最后会趋近于K，而且T越大，图像越平缓，T越小，图像越陡。