當前位置:編程學習大全網 - 編程語言 - 用MTALAB編程實現四階龍格庫塔解二元二階微分方程組

用MTALAB編程實現四階龍格庫塔解二元二階微分方程組

#include<stdlib.h>

#include<stdio.h>

/*n表示幾等分,n+1表示他輸出的個數*/

int RungeKutta(double y0,double a,double b,int n,double *x,double *y,int style,double (*function)(double,double))

{

double h=(b-a)/n,k1,k2,k3,k4;

int i;

// x=(double*)malloc((n+1)*sizeof(double));

// y=(double*)malloc((n+1)*sizeof(double));

x[0]=a;

y[0]=y0;

switch(style)

{

case 2:

for(i=0;i<n;i++)

{

x[i+1]=x[i]+h;

k1=function(x[i],y[i]);

k2=function(x[i]+h/2,y[i]+h*k1/2);

y[i+1]=y[i]+h*k2;

}

break;

case 3:

for(i=0;i<n;i++)

{

x[i+1]=x[i]+h;

k1=function(x[i],y[i]);

k2=function(x[i]+h/2,y[i]+h*k1/2);

k3=function(x[i]+h,y[i]-h*k1+2*h*k2);

y[i+1]=y[i]+h*(k1+4*k2+k3)/6;

}

break;

case 4:

for(i=0;i<n;i++)

{

x[i+1]=x[i]+h;

k1=function(x[i],y[i]);

k2=function(x[i]+h/2,y[i]+h*k1/2);

k3=function(x[i]+h/2,y[i]+h*k2/2);

k4=function(x[i]+h,y[i]+h*k3);

y[i+1]=y[i]+h*(k1+2*k2+2*k3+k4)/6;

}

break;

default:

return 0;

}

return 1;

}

double function(double x,double y)

{

return y-2*x/y;

}

//例子求y'=y-2*x/y(0<x<1);y0=1;

/*

int main()

{

double x[6],y[6];

printf("用二階龍格-庫塔方法\n");

RungeKutta(1,0,1,5,x,y,2,function);

for(int i=0;i<6;i++)

printf("x[%d]=%f,y[%d]=%f\n",i,x[i],i,y[i]);

printf("用三階龍格-庫塔方法\n");

RungeKutta(1,0,1,5,x,y,3,function);

for(i=0;i<6;i++)

printf("x[%d]=%f,y[%d]=%f\n",i,x[i],i,y[i]);

#include<stdlib.h>

#include<stdio.h>

/*n表示幾等分,n+1表示他輸出的個數*/

int RungeKutta(double y0,double a,double b,int n,double *x,double *y,int style,double (*function)(double,double))

{

double h=(b-a)/n,k1,k2,k3,k4;

int i;

// x=(double*)malloc((n+1)*sizeof(double));

// y=(double*)malloc((n+1)*sizeof(double));

x[0]=a;

y[0]=y0;

switch(style)

{

case 2:

for(i=0;i<n;i++)

{

x[i+1]=x[i]+h;

k1=function(x[i],y[i]);

k2=function(x[i]+h/2,y[i]+h*k1/2);

y[i+1]=y[i]+h*k2;

}

break;

case 3:

for(i=0;i<n;i++)

{

x[i+1]=x[i]+h;

k1=function(x[i],y[i]);

k2=function(x[i]+h/2,y[i]+h*k1/2);

k3=function(x[i]+h,y[i]-h*k1+2*h*k2);

y[i+1]=y[i]+h*(k1+4*k2+k3)/6;

}

break;

case 4:

for(i=0;i<n;i++)

{

x[i+1]=x[i]+h;

k1=function(x[i],y[i]);

k2=function(x[i]+h/2,y[i]+h*k1/2);

k3=function(x[i]+h/2,y[i]+h*k2/2);

k4=function(x[i]+h,y[i]+h*k3);

y[i+1]=y[i]+h*(k1+2*k2+2*k3+k4)/6;

}

break;

default:

return 0;

}

return 1;

}

printf("用四階龍格-庫塔方法\n");

RungeKutta(1,0,1,5,x,y,4,function);

for(i=0;i<6;i++)

printf("x[%d]=%f,y[%d]=%f\n",i,x[i],i,y[i]);

return 1;

  • 上一篇:CSS——知識點補充(四)元素的浮動屬性
  • 下一篇:我是學電子信息工程的,想向計算機網絡方向發展,該怎麽學,學什麽
  • copyright 2024編程學習大全網