龙格库塔算法解微分方程组 c语言
龙格库塔法解微分方程组
龙格库塔法解微分方程组引言微分方程组是数学中经常遇到的问题，在物理、工程和自然科学中都有广泛的应用。
为了求解微分方程组，我们需要利用数值方法来逼近解析解。
本文将介绍一种常用的数值方法——龙格库塔法（Runge-Kutta method），并探讨如何利用该方法来解微分方程组。
龙格库塔法概述龙格库塔法是一种迭代数值方法，用于求解常微分方程的初值问题。
它的主要思想是将微分方程的解进行离散化，将其转化为一系列的逼近值。
龙格库塔法的基本步骤如下：1.确定步长h和迭代次数n。
2.初始化初始条件，并假设第一个逼近值为y(xi)。
3.依次计算每个逼近值，直到得到y(xi+n*h)为止。
在每个迭代步骤中，龙格库塔法根据前一步的逼近值来计算下一步的逼近值。
该方法具有高阶精度和较好的稳定性，在实际应用中广泛使用。
单一微分方程的龙格库塔法首先，我们来看如何利用龙格库塔法来解一阶常微分方程。
以方程dy/dx = f(x, y)为例，其中f(x, y)为给定的函数。
步骤一：确定步长和迭代次数选择合适的步长h和迭代次数n来进行数值计算。
步长h决定了离散化的精度，而迭代次数n决定了逼近解的数目。
步骤二：初始化条件并计算逼近值设初始条件为y(x0) = y0，其中x0为起始点，y0为起始点处的函数值。
我们先通过欧拉法计算出y(x0 + h)的逼近值，然后再通过该逼近值来计算下一个逼近值。
逼近值的计算公式如下：k1 = h * f(x0, y0)k2 = h * f(x0 + h/2, y0 + k1/2)k3 = h * f(x0 + h/2, y0 + k2/2)k4 = h * f(x0 + h, y0 + k3)y(x0 + h) = y0 + 1/6 * (k1 + 2k2 + 2k3 + k4)步骤三：重复步骤二直到得到y(xi+n*h)依次利用上一步计算出的逼近值来计算下一个逼近值，直到得到y(xi+n*h)为止。
微分方程组的龙格库塔法对于一阶微分方程组的初值问题，我们可以将其转化为向量形式。
龙格库塔解答
实验题目3 四阶Runge-Kutta 方法摘要00(,)()dyf x y dxy x y⎧=⎪⎨⎪=⎩（6.1） 常微分方程初值问题的数值解法是求方程（6.1）的解在点列1(0,1,)n n n x x h n -=+= 上的近似值n y ，这里n h 是1n x -到n x 的步长，一般略去下标记为h 。
经典的R K -方法是一个四阶的方法，它的计算公式是：112341213243(22)6(,)(,)22(,)22(,)n n n n n n n n n n h y y K K K K K f x y h h K f x y K h h K f x y K K f x h y hK +⎧=++++⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩ （6.2）R K-方法的优点是：单步法、精度高，计算过程便于改变步长，缺点是计算量较大，每前进一步需要计算四次函数值f 。
在用龙格库塔方法时,要注意N 的选择要合适，N 太大，会使计算量加大，N 太小，h 较大，可能会使误差增大。
因此选择合适的N 很重要。
我们要在考虑精度的基础上，选择合适的N 。
在此，用c 语言实现了龙格库塔方法。
前言在此程序中把N 定义成宏常量，运行该程序，只需输入x0，y0，修改函数f ，再运行程序即可。
龙格-库塔法流程图程序设计流程#include<stdio.h>#define N 5double f(double x,double y);void main(){double a=0,b,aa,h,x,y,K1,K2,K3,K4;long n;printf("please input a,b,aa\n");scanf("%lf,%lf,%lf",&a,&b,&aa);h=(double)(b-a)/N;for(n=1;n<=N;n++){K1=h*f(a,aa);K2=h*f(a+h/2,aa+K1/2);K3=h*f(a+h/2,aa+K2/2);K4=h*f(a+h,aa+K3);x=a+h;y=aa+(double)1/6*(K1+2*K2+2*K3+K4);printf("x1=%lf,y1=%lf\n",x,y);a=x;aa=y;}}double f(double x,double y){double f;f=x+y;return f;}龙格库塔思考题解答:问题1解答:数值解和解析解不一样,因为数值解是对给定的x值求出y值,而没有给出y与x的函数关系.相反解析解就给出了y与x的函数关系。
龙格—库塔 C++
实验四 龙格-库塔法【实验内容】1. 用标准四阶龙格－库塔方法设计算法;2. 编程解微分方程初值问题.3.在课后习题中选择一个题目编程计算。
交回实验报告与计算结果。
【实验方法或步骤】四阶龙格－库塔方法的计算步骤求解'0()(,)()()y x f x y a x b y a y ⎧=≤≤⎨=⎩ 对上述给定的(,)f x y ，用四阶龙格-库塔法求解常微分方程初值问题112341213243(22)6(,)11(,)2211(,)22(,)n n n n n n n n n n h y y k k k k k hf x y k hf x h y k k hf x h y k k hf x h y k +⎧=++++⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩[实验程序]（自编程序 .cpp 文件）#include<stdio.h>double f(double a, double b){double c;c = a + b;return c;}int main(void)inti;double a, b, m, n, h;double k1, k2, k3, k4;printf("请输入端点值a和b：\n");scanf("%lf %lf", &a, &b);printf("请输入整数n：\n");scanf("%lf", &n);printf("请输入初始值：\n");scanf("%lf", &m);h = (b-a)/n;printf("在区间（%3.1lf %3.1lf）上步长为%3.1lf初值为: \n", a, b, h);printf("X0 = %lf, Y0 = %lf\n", a ,m);for (i=0; i<n; i++){k1 = f(a, m);k2 = f(a+h/2, m+h*k1/2);k3 = f(a+h/2, m+h*k2/2);k4 = f(a+h, m+h*k3);m = m + (k1+2*k2+2*k3+k4)*h/6;a = a + h;printf("X%d = %lf, Y%d = %lf\n",i+1, a ,i+1, m);}return 0;【实验结果】（课后4题）。
C语言 微分方程组龙格-库塔法
/* By wende Isaac微分方程组龙格－库塔法*/#include<stdio.h>#include<math.h>double Ft(int j, double x ,double y[]);void rkutta4g(double x0, double y0[], double ym[], int M,double h);int main(){int N=5000,i,j,ipm=0,ip=100,M=50;double H=0.05;double x[100],y[100][10];double x0,y0[M],ym[M];x[0]=0.0;y[0][0]=1.0;for(i=1;i<M;i++){y[0][i]=0;}for(i=1;i<=N;i++){x[i]=x[i-1]+H;ipm++;x0=x[i-1];for(j=0;j<M;j++){y0[j]=y[i-1][j];}double ym[M];rkutta4g(x0,y0,ym,M,H);for(j=0;j<M;j++){y[i][j]=ym[j];// printf("ym[%d]= %f \n",j,ym[j]);}if(ipm>=ip){printf(" T: %.f(sec) ",x[i]);for (j=0;j<M-1;j++){printf(" C%d= %0.5f(m) ",j+1,y[i][j]);}printf(" C%d= %f(m) \n",M,y[i][j]);ipm=0;}}return 0;}void rkutta4g(double x0, double y0[], double ym[], int M,double h) {double s[4]={0,0.5,0.5,1},k[M][5],c[4]={1.0/6,1.0/3,1.0/3,1.0/6};int i,j,l;for(i=0;i<M;i++){k[i][0]=0;}for(j=0;j<M;j++){ym[j]=y0[j];// printf("ym[%d]=%f\n",j,ym[j]);}double xt,yt[5];for(j=0;j<4;j++){xt=x0+s[j]*h; //printf("xt=%f \n",xt);for(i=0;i<M;i++){yt[i]=y0[i]+s[j]*k[i][j];// printf("yt[i]=%f \n",yt[i]);}for(i=0;i<M;i++){k[i][j+1]=h*Ft(i,xt,yt);// printf("k[%d][%d] =%f \n",i,j+1,k[i][j+1]);}for(i=0;i<M;i++){ym[i]=ym[i]+c[j]*k[i][j+1];}}}double Ft(int j, double x ,double y[]){double k0,k1,k2,k3,ft;switch(j){case 0:ft=导数方程1;break;case 1:ft=导数方程2;break;case 2:ft=导数方程13;break;case 3:ft=导数方程4;break;case 4:ft=导数方程5;break;default :printf("error operation \n");}// printf("ft= %f \n", ft);return ft;}。
龙格库塔法解常微分方程
龙格库塔法解常微分方程龙格库塔法是一种常用的数值解决常微分方程（ODE）的数值计算方法。
它于1960年末由古斯塔夫·龙格库塔（Gustav Runge-Kutta）提出，并多次得到不同的改进和拓展，成为经典的数值解决ODE的一种方法。
龙格库塔法的基本思想是，将ODE的微分方程化为一组非线性的简单方程，通过求解这些简单方程来近似解决ODE。
龙格库塔法比较通用，可以解决一阶和多阶的常微分方程，其代表性的多阶龙格库塔方程有隐式四阶龙格库塔方程和显式四阶龙格库塔方程。
设微分方程组如下：$$\frac{ d y}{ d x }=f(x,y(x))$$由此可以推导出隐式四阶龙格库塔方程（前向差分形式）：$$y_{n+1}=y_{n}+\frac{h}{6}(k_{1}+2k_{2}+2k_{3}+k_{4}) $$其中$h$是步长，$k_{1}=f(x_{n},y_{n})$，$k_{2}=f(x_{n}+\frac{h}{2},y_{n}+\frac{h}{2}k_{1})$，$k_{3}=f(x_{n}+\frac{h}{2},y_{n}+\frac{h}{2}k_{2})$，$k_{4}=f(x_{n}+h,y_{n}+hk_{3})$。
龙格库塔法有一个重要特点，就是具有“四步”方案，即在一个步长下只计算四次，完成整个迭代计算过程。
这使得龙格库塔法可以获得高效、稳定可靠的数值解。
补充来说，不同的函数问题有不同的龙格库塔方程，比如一阶非线性龙格库塔方程、二阶线性龙格库塔方程等。
总的说来，龙格库塔法是一种用来数值解常微分方程的有效、可靠的方法，它可以有效解决一阶、多阶和非线性的微分方程。
以分析性形式构建的方程组都可以得到精确的数值解，而且计算量也比较少，速度较快，是一种常用的数值求解ODE的方法。
因此，龙格库塔法作为一种能够有效解决常微分方程的数值计算方法，在理论和实践中均受到广泛的应用。
相比传统的求解方法，它具有更高的计算效率，能够更快速、更准确地给出ODE的数值解。
四阶龙格库塔法求解微分方程组
四阶龙格库塔法求解微分方程组微分方程是数学中的一个重要分支，它描述了自然界中许多现象的发展规律。
在实际应用中，我们经常需要求解微分方程组，以得到系统的演化规律和性质。
本文将介绍一种常用的求解微分方程组的数值方法——四阶龙格库塔法。
一、微分方程组的数值解法微分方程组是形如下式的方程集合：begin{cases} frac{dy_1}{dx}=f_1(x,y_1,y_2,cdots,y_n) frac{dy_2}{dx}=f_2(x,y_1,y_2,cdots,y_n) cdotsfrac{dy_n}{dx}=f_n(x,y_1,y_2,cdots,y_n) end{cases} 其中，$y_1,y_2,cdots,y_n$是关于自变量$x$的未知函数，$f_1,f_2,cdots,f_n$是已知的函数。
求解微分方程组的数值方法主要有以下两种：1.欧拉法欧拉法是最简单的数值方法之一，它的基本思想是将微分方程组中的导数用差分代替，从而得到一个递推公式。
具体而言，欧拉法的递推公式为：$$y_{i+1}=y_i+hf(x_i,y_i)$$其中，$h$是步长，$x_i$和$y_i$分别是第$i$个点的自变量和因变量的值。
欧拉法的精度较低，容易出现数值误差，但是它的计算速度快，适用于一些简单的微分方程组。
2.龙格库塔法龙格库塔法是一种常用的高精度数值方法，它的基本思想是将微分方程组中的导数用一系列加权的差分代替，从而得到一个递推公式。
其中，四阶龙格库塔法是最为常用的一种。
具体而言，四阶龙格库塔法的递推公式为：begin{aligned} k_1&=hf(x_i,y_i)k_2&=hf(x_i+frac{h}{2},y_i+frac{k_1}{2})k_3&=hf(x_i+frac{h}{2},y_i+frac{k_2}{2})k_4&=hf(x_i+h,y_i+k_3)y_{i+1}&=y_i+frac{k_1+2k_2+2k_3+k_4}{6} end{aligned} 其中，$k_1,k_2,k_3,k_4$是加权的差分，$h$是步长，$x_i$和$y_i$分别是第$i$个点的自变量和因变量的值。
c求解四阶龙格库塔法
某一地区的病菌传染，三种人员的人数的状态方程，即可能受传染的人数x 1,已被传染的病的人数x2，及已治愈的人数x3，并假设他们是因接触而传染。
设：α是x1中单位时间内的传染系数，β是x2中单位时间内被治愈的比例系数，则得以下状态方程：f`(x1)=-α*x1*x2;f`(x2)=α*x1*x2-β*x2;f`(x3)=β*x2;其中f`(x)是加速度本题中，α=0.01，β=0.027初始条件：t=0, x1=620人, x2=10人, x3=70人选步长：h=0.025, 用四阶RK法模拟。
-------------程序清单-------------------#include <stdio.h>#include <conio.h>#define A 0.01#define B 0.027void myFun(double *x,double *f){f[0]=-A*x[0]*x[1];f[1]=A*x[0]*x[1]-B*x[1];f[2]=B*x[1] ;}int rht4(double *rs,double h,double x[],int n,int Num){int j,i,ii;double ps[4]={0,0.5,0.5,1};double bs[4]={1.0/6,1.0/3,1.0/3,1.0/6};double K[3];double y[3];double f[3];for(i=0;i<n;i++){y[i]=x[i];rs[i]=y[i];}for(j=1;j<Num;j++){for(ii=0;ii<4;ii++){ /*龙格-库塔*/for(i=0;i<n;i++){x[i]=rs[(j-1)*n+i]+ps[ii]*f[i]*h; /*特别注意rs[(j-1)*n+i] 而不是x[i] 不是前次的*/}myFun(x,f);/*算每次的加速度*/for(i=0;i<n;i++){y[i]+=f[i]*h*bs[ii];}} /*end 龙格-库塔*/for(i=0;i<n;i++){x[i]=y[i];rs[j*n+i]=y[i];}}}int main(){double t=0,x[3]={620.0,10.0,70.0},h=0.025,*rh;int Num,i,j;double f[3];printf("Input Num: ");scanf("%d",&Num);rh=(double *)malloc(Num*3*sizeof(double));rht4(rh,h,x,3,Num) ;printf("x1 x2 x3-----------------------------------------------------"); for(i=0;i<Num;i++){printf(" ");for(j=0;j<3;j++)printf("%f ",rh[i*3+j]);}getch();return 0;}－－－－－－－－－－end－－－－－－－－－－－－－－－。
龙格库塔算法解微分方程组c语言
/*＊＊*＊*＊＊＊＊*****＊＊*＊＊＊**＊＊***＊＊****＊*＊＊*************＊＊*＊＊*＊＊*＊**＊*＊*＊*＊＊*＊*＊This program is to solve the initial value problem of following systemof differential equations：dx/dt=x+2*y，x（0）=0，dy/dt=2＊x+y,y(0)=2，x(0.2) and y(0。
2） are to be calculated＊*＊＊**＊*＊*＊*＊**＊＊**＊＊＊****＊＊＊*＊*＊*＊*****＊＊***＊＊＊＊＊*＊＊**＊＊**＊*＊＊＊*＊***＊＊＊＊＊*＊/#include<stdio。
h>＃include<math。
h〉#define steplength 0。
1 //步?长¡è可¨¦根¨´据Y需¨¨要°a调Ì¡Â整？;＃define FuncNumber 2 //FuncNumber为a微¡é分¤?方¤?程¨¬的Ì?数ºy目？;void main(){float x[200]，Yn［20][200］，reachpoint；int i；x[0］=0；Yn[0][0］=0；Yn［1］[0］=2； //初？值¦Ì条¬？件t；reachpoint=0。
2; //所¨´求¨®点Ì？可¨¦根¨´据Y需¨¨要°a调Ì¡Â整？;void rightfunctions(float x ,float *Auxiliary，float *Func）;void Runge_Kutta(float *x，float reachpoint， float(＊Yn）[200］)；Runge_Kutta（x ，reachpoint， Yn);printf("x ")；for(i=0;i<=（reachpoint—x[0])/steplength；i++)printf("%f "，x[i]);printf("\nY1 ”)；for（i=0；i〈=（reachpoint—x［0])/steplength;i++)printf("％f ”,Yn［0]［i］）；printf(”\nY2 ”）；for(i=0;i〈=（reachpoint—x［0］）/steplength;i++）printf("%f ”，Yn[1]［i］)；getchar（）；}void rightfunctions（float x ,float *Auxiliary,float *Func）//当Ì¡À右®¨°方¤?程¨¬改?变À？时º¡À,ê？需¨¨要°a改?变À？；{Func［0］=Auxiliary[0]+2*Auxiliary[1］;Func［1］=2*Auxiliary［0］+Auxiliary[1]；}void Runge_Kutta（float *x,float reachpoint, float（＊Yn)［200］)｛int i，j;float Func［FuncNumber］,K［FuncNumber]［4]，Auxiliary[FuncNumber］；for(i=0；i<=(reachpoint—x[0])/steplength;i++)｛for（j=0;j<FuncNumber；j++)Auxiliary[j]=＊(Yn［j］+i)；rightfunctions（x[i］，Auxiliary，Func）;for(j=0；j<FuncNumber;j++）｛K[j][0]=Func[j];Auxiliary[j］=*（Yn[j］+i)+0.5＊steplength＊K［j］［0］；}rightfunctions(x[i]，Auxiliary,Func）；for（j=0;j<FuncNumber;j++)｛K［j］[1］=Func[j］;Auxiliary［j]=＊（Yn[j]+i)+0.5＊steplength*K［j］［1］；}rightfunctions（x[i］,Auxiliary，Func）;for（j=0;j〈FuncNumber;j++){K［j]［2］=Func［j];Auxiliary[j］=＊(Yn[j]+i)+steplength＊K［j]［2］；｝rightfunctions（x[i］,Auxiliary,Func)；for(j=0；j〈FuncNumber;j++)K[j]［3]=Func［j］；for（j=0；j〈FuncNumber;j++)Yn[j]［i+1］=Yn[j］[i］+（K[j][0］+2*K[j］[1］+2*K[j]［2]+K［j][3］)*steplength/6.0；x[i+1]=x[i］+steplength；｝}。
龙格库塔法的c
龙格库塔法的c++编程CODE:#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]);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;龙格库塔求解微分方程数值解工程中很多的地方用到龙格库塔求解微分方程的数值解，龙格库塔是很重要的一种方法，尤其是四阶的，精确度相当的高。
使用C语言解常微分方程 C ODE
解常微分方程姓名：Vincent年级：2010，学号：1033****，组号：5（小组），4（大组）1. 数值方法:我们的实验目标是解常微分方程，其中包括几类问题。
一阶常微分初值问题，高阶常微分初值问题，常微分方程组初值问题，二阶常微分方程边值问题，二阶线性常微分方程边值问题。
对待上面的几类问题，我们分别使用不同的方法。
∙ 初值问题使用 龙格-库塔 来处理 ∙ 边值问题用打靶法来处理 ∙ 线性边值问题有限差分法初值问题我们分别使用∙ 二阶 龙格-库塔 方法 ∙ 4阶 龙格-库塔 方法 来处理一阶常微分方程。
理论如下：对于这样一个方程'()(,)y t f t y =当h 很小时，我们用泰勒展开，12111111(,)(,)(,)k k k k i i k k j j i ij k hf t y k hf t h y k k hf a b a b t h y h k -===++=++∑当我们选择正确的参数 a[i][j],b[i][j]之后，就可以近似认为这就是泰勒展开。
对于二阶，我们有：()01()()y t h y t h Af Bf+=++其中010(,)(,)P f f t y f f t y h Qhf ==++经过前人的计算机经验，我们有，选择 A=1/2,B=1/2,则 P=1,Q=1，于是又如下形式，这也使休恩方法的表达式。
所以我们称其为 龙格（库塔）休恩方法()()()(,)(,(,))2hy t h y t f t y f t h y hf t y +=++++对于4阶龙格方法，我们有类似的想法， 我们使用前人经验的出的系数，有如下公式112341213243(22),6(,),(,),22(,),22(,).n n n n n n n n n n h y y k k k k k f t y h h k f t y k h h k f t y k k f t h y hk +⎧=++++⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩对于高阶微分方程及微分方程组 我们用∙ 4阶龙格-库塔方法来解对于一个如下的微分方程组'111'1(,,,),(,,,).n nn n y f t y y y f t y y ⎧=⎪⎨⎪=⎩ 101,00,0(),()n n y t y y t y=⎧⎪⎨⎪=⎩ 我们可以认为是一个一阶向量微分方程，所以可以用龙格-库塔方法的向量形式解。
经典四阶龙格库塔法解一阶微分方程组
数值计算课程设计1. 经典四阶龙格库塔法解一阶微分方程组1.1 运用四阶龙格库塔法解一阶微分方程组算法分析x k 1 x k h(f 1 2 f 2 2 f 3 f 4)6hy k 1 y k h 6(g 1 2g 2 2g 3 g 4 )t k 1 t k h经过循环计算由 t 0,x 0, y 0推得 t 1,x 1,y1 t 2,x 2,y 2⋯⋯每个龙格 -库塔方法都是由一个合适的泰勒方法推导而来，使得其最终全局 误差为 O h N, 一种折中方法是每次进行若干次函数求值，从而省去高阶导数计 算。
4阶龙格-库塔方法 (RK4)是最常用的，它适用于一般的应用，因为它非常精 准，稳定，且易于编程。
f 1 f (t k , x k , y k ) ，hf 2 f (t k h 2,x khh2f 1,y k 2g 1)hf 3 f (t k h 2,x khf 2,y k 2g 2)f 4 f (t k h ,x k hf 3,y k hg 3)g 1 g(t k ,x k , y k ) h g 2g(t k 2,x khh2f 1,y k 2g 1)hg 3 g(t k 2, x khf 2h2,y k 2g 2) g 4 g(t k h,x k hf 3, y k hg 3)1-1)1-2)1-3)1-4)1-5)1-6)1-7)1-8)1-9)1-10 )经典四阶龙格库塔法解一阶微分方程组1.2 经典四阶龙格库塔法解一阶微分方程流程图1.3 经典四阶龙格库塔法解一阶微分方程程序代码：#include <iostream> #include <iomanip> using namespace std;void RK4( double (*f)(double t,double x, double y),double (*g)(double t,double x, double y) ,double initial[3], double resu[3],double h){double f1,f2,f3,f4,g1,g2,g3,g4,t0,x0,y0,x1,y1; t0=initial[0];x0=initial[1];y0=initial[2]; f1=f(t0,x0,y0); g1=g(t0,x0,y0); f2=f(t0+h/2, x0+h*f1/2,y0+h*g1/2);g2=g(t0+h/2,x0+h*f1/2,y0+h*g1/2); f3=f(t0+h/2,x0+h*f2/2,y0+h*g2/2);g3=g(t0+h/2,x0+h*f2/2,y0+h*g2/2);f4=f(t0+h, x0+h*f3,y0+h*g3); g4=g(t0+h, x0+h*f3,y0+h*g3);图 1-1 经典四阶龙格库塔法解一阶微分方程流程图数值计算课程设计x1=x0+h*(f1+2*f2+2*f3+f4)/6; y1=y0+h*(g1+2*g2+2*g3+g4)/6;resu[0]=t0+h;resu[1]=x1;resu[2]=y1;}int main(){double f(double t,double x, double y);double g(double t,double x, double y);double initial[3],resu[3];double a,b,H;double t,step;int i;cout<<" 输入所求微分方程组的初值t0,x0,y0:"; cin>>initial[0]>>initial[1]>>initial[2]; cout<<" 输入所求微分方程组的微分区间[a,b]:"; cin>>a>>b;cout<<" 输入所求微分方程组所分解子区间的个数step:"; cin>>step;cout<<setiosflags(ios::right)<<setiosflags(ios::fixed)<<setprecision( 10);H=(b-a)/step;cout<< initial[0]<<setw(18)<<initial[1]<<setw(18)<<initial[2]<<endl;for(i=0;i<step;i++){ RK4( f,g ,initial, resu,H);cout<<resu[0]<<setw(20)<<resu[1]<<setw(20)<<resu[2]<<endl;initial[0]=resu[0];initial[1]=resu[1];initial[2]=resu[2];}return(0);}double f(double t,double x, double y){double dx;经典四阶龙格库塔法解一阶微分方程组dx=x+2*y; return(dx);}double g(double t,double x, double y) {double dy;dy=3*x+2*y; return(dy);}1.4 经典四阶龙格库塔法解一阶微分方程程序调试结果图示：图 1-2 经典四阶龙格库塔法解一阶微分方程算法程序调试图2. 高斯列主元法解线性方程组2.1 高斯列主元法解线性方程组算法分析使用伪代码编写高斯消元过程：for k=1 to n-1 dofor i=k+1 to n l<=a(i,k)/a(k,k) for j=k to n do a(i,j)<=a(i,j)-l*a(k,j) end %end of for j b(i)<=b(i)-l*b(k) end %end offor i end %end of for k最后得到 A ，b 可以构成上三角线性方程组 接着使用回代法求解上三角线性方程组应用所编写程序计算所给例题：其中初值为求解区间为 。
四阶龙格库塔法求解微分方程组
四阶龙格库塔法求解微分方程组四阶龙格-库塔法（RK4）是一种常用的数值求解微分方程组的方法。
它通过将微分方程转化为差分方程，通过逐步迭代来求得方程的近似解。
首先，考虑一个一阶微分方程组：dx/dt = f(x,t)其中，x是一个向量，f是一个向量函数，t是时间。
我们的目标是求解x(t)。
RK4法的基本思想是，在每个步骤中，通过预测下一个时间点的解并进行修正来逼近实际的解。
具体步骤如下：1.给定初始条件x(t0)和步长h，计算t=t0+h。
2. 计算k1 = hf(x(t0),t0)，即利用初始点(x(t0),t0)计算出的斜率。
3. 计算k2 = hf(x(t0)+k1/2,t0+h/2)，即利用中间点(x(t0)+k1/2,t0+h/2)计算出的斜率。
4. 计算k3 = hf(x(t0)+k2/2,t0+h/2)，同理得到另一个中间点的斜率。
5. 计算k4 = hf(x(t0)+k3,t0+h)，即利用最后一个中间点(x(t0)+k3,t0+h)计算的斜率。
6.根据k1，k2，k3和k4计算下一个点的值：x(t1)=x(t0)+(k1+2k2+2k3+k4)/67.重复步骤2到6，直到达到所需的终止时间。
接下来，让我们通过一个具体的例子来说明RK4方法的应用。
假设我们要求解如下的一阶微分方程组：dx/dt = f(x,y,z)dy/dt = g(x,y,z)dz/dt = h(x,y,z)其中f，g和h是已知的向量函数，我们需要找到x(t)，y(t)和z(t)。
首先，我们需要将该微分方程组转化为差分方程组。
我们将离散时间点t0，t1，t2，...，tn代入微分方程，得到：x(t0+h)=x(t0)+h*f(x(t0),y(t0),z(t0))y(t0+h)=y(t0)+h*g(x(t0),y(t0),z(t0))z(t0+h)=z(t0)+h*h(x(t0),y(t0),z(t0))然后，我们通过逐步迭代的方式来求解方程组。
经典四阶龙格库塔法解一阶微分方程组
1.经典四阶龙格库塔法解一阶微分方程组1.1运用四阶龙格库塔法解一阶微分方程组算法分析),,(1k k k y x t f f =， )2,2,2(112g hy f h x h t f f k k k +++=)2,2,2(223g hy f h x h t f f k k k +++=),,(334hg y hf x h t f f k k k +++=),,(1k k k y x t g g =)2,2,2(112g hy f h x h t g g k k k +++=)2,2,2(223g hy f h x h t g g k k k +++=),,(334hg y hf x h t g g k k k +++=)22(6)22(64321143211g g g g hy y f f f f hx x k k k k ++++=++++=++1k k t t h +=+经过循环计算由 推得 ……每个龙格-库塔方法都是由一个合适的泰勒方法推导而来，使得其最终全局误差为()N O h ,一种折中方法是每次进行若干次函数求值，从而省去高阶导数计算。
4阶龙格-库塔方法(RK4)是最常用的，它适用于一般的应用，因为它非常精准，稳定，且易于编程。
000,,t x y ()()111222,,,,t x y t x y （1-1）（1-2）（1-3）（1-4）（1-5）（1-6）（1-7）（1-8）（1-9）（1-10）1.2经典四阶龙格库塔法解一阶微分方程流程图图1-1 经典四阶龙格库塔法解一阶微分方程流程图1.3经典四阶龙格库塔法解一阶微分方程程序代码：#include <iostream>#include <iomanip>using namespace std;void RK4( double (*f)(double t,double x, double y),double (*g)(double t,double x, double y) ,double initial[3], double resu[3],double h) {double f1,f2,f3,f4,g1,g2,g3,g4,t0,x0,y0,x1,y1;t0=initial[0];x0=initial[1];y0=initial[2];f1=f(t0,x0,y0); g1=g(t0,x0,y0);f2=f(t0+h/2, x0+h*f1/2,y0+h*g1/2); g2=g(t0+h/2, x0+h*f1/2,y0+h*g1/2);f3=f(t0+h/2, x0+h*f2/2,y0+h*g2/2); g3=g(t0+h/2,x0+h*f2/2,y0+h*g2/2);f4=f(t0+h, x0+h*f3,y0+h*g3); g4=g(t0+h, x0+h*f3,y0+h*g3); x1=x0+h*(f1+2*f2+2*f3+f4)/6; y1=y0+h*(g1+2*g2+2*g3+g4)/6; resu[0]=t0+h;resu[1]=x1;resu[2]=y1;}int main(){double f(double t,double x, double y);double g(double t,double x, double y);double initial[3],resu[3];double a,b,H;double t,step;int i;cout<<"输入所求微分方程组的初值t0,x0,y0:";cin>>initial[0]>>initial[1]>>initial[2];cout<<"输入所求微分方程组的微分区间[a,b]:";cin>>a>>b;cout<<"输入所求微分方程组所分解子区间的个数step:";cin>>step;cout<<setiosflags(ios::right)<<setiosflags(ios::fixed)<<setprecision( 10);H=(b-a)/step;cout<< initial[0]<<setw(18)<<initial[1]<<setw(18)<<initial[2]<<endl; for(i=0;i<step;i++){ RK4( f,g ,initial, resu,H);cout<<resu[0]<<setw(20)<<resu[1]<<setw(20)<<resu[2]<<endl;initial[0]=resu[0];initial[1]=resu[1];initial[2]=resu[2];}return(0);}double f(double t,double x, double y){double dx;dx=x+2*y; return(dx); }double g(double t,double x, double y) {double dy; dy=3*x+2*y; return(dy); }1.4经典四阶龙格库塔法解一阶微分方程程序调试结果图示：应用所编写程序计算所给例题：其中初值为求解区间为[0,0.2]。
龙格-库塔求解微分方程
MATLAB 在数值分析中的应用摘 要：数值分析是研究数学问题数值解及其理论的一个数学分支,涉及面很广.自计算机在其中应用以后,使其应用价值更为广泛,解决实际问题更为有效.以MA TLAB 为平台的数值分析方法与图形可视化编程在解决一些复杂问题时,大大降低了难度,减少了工作量.本文以经典龙格-库塔算法为核心,介绍用数值方法解常微分方程(组),分析自制系统的稳定性问题.如范德波尔方程.以及混沌现象中奇怪吸引子的模拟实现关键词： MA TLAB; 经典龙格-库塔算法; 解常微分方程(组); 范德波尔方程; 奇怪吸引子;1 MATLAB 平台1.1 矩阵MATLAB 是矩阵实验室（Matrix Laboratory ）的简称, 美国MathWorks 公司出品的商业数学软件.以下介绍在处理方程组时的特殊技法: A=[a b c; d e f; g h i]得到矩阵a b c d e f g h i也可直接输入方程组:function y=equition(t,x)y=[q(2)+q(1)/sqrt(q(1)^2+q(2)^2)*(1-q(1)^2-q(2)^2); % 以q 向量保存方程组;-q(1)+q(2)/sqrt(q(1)^2+q(2)^2)*(1-q(1)^2-q(2)^2)] % q(1)代表应变量x, q(2) 代表应变量y得到1.2 语言MATLAB 编程语言与 C 语言 类似, 格式上略有不同. 里不做详述, 程序中会做必要的标注.但在处理矩阵式十分方便, C++是以数组的形式存放变量,对某个元素都要准确指出其位置.22222222()]()]x y x y x y y x x y x y ⎧=+-+⎪+⎪⎨⎪=-+-+⎪+⎩C++是以数组的形式存放变量,对某个元素都要准确指出其位置,处理循环时算法复杂,即使使用指针也很容易出错.2经典龙格-库塔算法公式(15张)龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。
四阶龙格库塔法解一阶二元微分方程
fp=fopen("sjy.txt","w");
fp1=fopen("sjykxy.txt","w");
fprintf(fp1,"k\tx1\tx2\tx3\ty1\ty2\ty3\n");
if(fp==NULL||fp1==NULL){
printf("Failed to open file.\n");
xi=x[i]+dx;
yi=y[i]+dy;
return (xi-b*yi+a)/c;
}
void main(){
double Kx[3][4],Ky[3][4],x[3]={1,2,3},y[3]={2,3,4},xavg,k=0;//定义x,y的初值；
double z[3]={0};
int i,j,m,n,S;
fprintf(fp,"\n");
//取第S步，即时间为S*h的x1，x2，x3，y1，y2，y3随k值的变化；
while(j==S)
{
for(n=0;n<3;++n)
{
fprintf(fp1,"\t%.3lf",x[n]);
}
for(n=0;n<3;++n)
{
fprintf(fp1,"\t%.3lf",y[n]);
int j;
double xi,yi;
xi=x[i]+dx;
yi=y[i]+dy;
return c*(xi-pow(xi,3)/3+yi)+k*(xavg-xi)+c*z[i];
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
/***************************************************************************
This program is to solve the initial value problem of following system
of differential equati ons:
dx/dt=x+2*y,x(0)=0,
dy/dt=2*x+y,y(0)=2,
x(0.2) and y(0.2) are to be calculated
****************************************************************************/
#include<stdi o.h>
#include<math.h>
#define steplength 0.1 //步?长¡è可¨¦根¨´据Y需¨¨要°a调Ì¡Â整?;
#define FuncNumber 2 //FuncNumber为a微¡é分¤?方¤?程¨¬的Ì?数ºy目?;
void main()
{
float x[200],Yn[20][200],reachpoint;i nt i;
x[0]=0;Yn[0][0]=0;Yn[1][0]=2; //初?值¦Ì条¬?件t;
reachpoint=0.2; //所¨´求¨®点Ì?可¨¦根¨´据Y需¨¨要°a调Ì¡Â整?;
void rightfunctions(float x ,float *Auxiliary,float *Func);
void R unge_Kutta(float *x,float reachpoint, float(*Yn)[200]);
Runge_Kutta(x ,reachpoi nt, Yn);
printf("x ");
for(i=0;i<=(reachpoint-x[0])/steplength;i++)
printf("%f ",x[i]);
printf("\nY1 ");
for(i=0;i<=(reachpoint-x[0])/steplength;i++)
printf("%f ",Yn[0][i]);
printf("\nY2 ");
for(i=0;i<=(reachpoint-x[0])/steplength;i++)
printf("%f ",Yn[1][i]);
getchar();
}
void rightfunctions(float x ,float *Auxiliary,float *Func)//当Ì¡À右®¨°方¤?程¨¬改?变À?时º¡À，ê?需¨¨要°a改?变À?;
{
Func[0]=Auxiliary[0]+2*Auxiliary[1];
Func[1]=2*Auxiliary[0]+Auxiliary[1];
}
void R unge_Kutta(float *x,float reachpoint, float(*Yn)[200])
{
int i,j;
float Func[FuncNumber],K[FuncNumber][4],Auxiliary[FuncNumber];
for(i=0;i<=(reachpoint-x[0])/steplength;i++)
{
for(j=0;j<FuncNumber;j++)
Auxiliary[j]=*(Yn[j]+i);
rightfunctions(x[i],Auxiliary,Func);
for(j=0;j<FuncNumber;j++)
{
K[j][0]=Func[j];
Auxiliary[j]=*(Yn[j]+i)+0.5*steplength*K[j][0];
}
rightfunctions(x[i],Auxiliary,Func);
for(j=0;j<FuncNumber;j++)
{
K[j][1]=Func[j];
Auxiliary[j]=*(Yn[j]+i)+0.5*steplength*K[j][1];
}
rightfunctions(x[i],Auxiliary,Func);
for(j=0;j<FuncNumber;j++)
{
K[j][2]=Func[j];
Auxiliary[j]=*(Yn[j]+i)+steplength*K[j][2];
}
rightfunctions(x[i],Auxiliary,Func);
for(j=0;j<FuncNumber;j++)
K[j][3]=Func[j];
for(j=0;j<FuncNumber;j++)
Yn[j][i+1]=Yn[j][i]+(K[j][0]+2*K[j][1]+2*K[j][2]+K[j][3])*steplength/6.0;
x[i+1]=x[i]+steplength;
}
}。