二维微分方程组 龙格库塔 数值解

#include <iostream>
#include "Possion.h"
#include <vector>
using namespace std;

typedef float (*func)(const float&, const float&, const float&, const float&, const float&, const float&, const float&, const int&, const int&);

void RungeKutta_2D(const float& t_init,const float& t_interval, const float& t_final, const float& x_init, const float& y_init, vector<float>& vec_t, vector<float>& vec_x, vector<float>& vec_y, func funcp1, func funcp2,
				   const float& lambda, const float& aerfa, const float& beita,const float& sourcespeed, const int& filesize, const int& leaverate){

	// curr;
	float t_cur = t_init;
	float x_cur = x_init;
	float y_cur = y_init;

	// next;
	float x_next, y_next, t_next;
	x_next = x_init;
	y_next = y_init;
	t_next = t_init;
	// container init;
	vec_t.push_back(t_init);
	vec_x.push_back(x_init);
	vec_y.push_back(y_init);

	// define the 4 value;
	float k1x, k2x, k3x, k4x;
	float k1y, k2y, k3y, k4y;

	// loop to the final;
	while(1){
		if(t_next <= t_final){
			//const float &lambda, const float& aerfa, const float& beita, const float& sourcespeed, const float& t, const float& x, const float& y, const int& filesize
			// cal k1,k2,k3,k4;
			k1x = funcp1(lambda, aerfa, beita, sourcespeed, t_cur, x_cur, y_cur, filesize, leaverate);
			k1y = funcp2(lambda, aerfa, beita, sourcespeed, t_cur, x_cur, y_cur, filesize, leaverate);
			k2x = funcp1(lambda, aerfa, beita, sourcespeed, t_cur + (t_interval / 2), x_cur + (t_interval / 2) * k1x, y_cur + (t_interval / 2) * k1y, filesize, leaverate);
			k2y = funcp2(lambda, aerfa, beita, sourcespeed, t_cur + (t_interval / 2), x_cur + (t_interval / 2) * k1x, y_cur + (t_interval / 2) * k1y, filesize, leaverate);
			k3x = funcp1(lambda, aerfa, beita, sourcespeed, t_cur + (t_interval / 2), x_cur + (t_interval / 2) * k2x, y_cur + (t_interval / 2) * k2y, filesize, leaverate);
			k3y = funcp2(lambda, aerfa, beita, sourcespeed, t_cur + (t_interval / 2), x_cur + (t_interval / 2) * k2x, y_cur + (t_interval / 2) * k2y, filesize, leaverate);
			k4x = funcp1(lambda, aerfa, beita, sourcespeed, t_cur + t_interval, x_cur + t_interval * k3x, y_cur + t_interval * k3y, filesize, leaverate);
			k4y = funcp2(lambda, aerfa, beita, sourcespeed, t_cur + t_interval, x_cur + t_interval * k3x, y_cur + t_interval * k3y, filesize, leaverate);
			// next;
			x_next = x_cur + (t_interval / 6) * (k1x + 2 * k2x + 2 * k3x + k4x);
			y_next = x_cur + (t_interval / 6) * (k1y + 2 * k2y + 2 * k3y + k4y);
			t_next = t_cur + t_interval;
			// append to result;
			vec_t.push_back(t_next);
			vec_x.push_back(x_next);
			vec_y.push_back(y_next);
			// keep data;
			t_cur = t_next;
			x_cur = x_next;
			y_cur = y_next;
		}else
			break;
	}
	return;
}

龙格库塔方法广泛应用在微分积分求数值解的领域，其精确度非常之高，理论相对复杂。

简单的说，是应用了级数的性质以及斜率的概念。上述代码是二维龙格库塔方法。采用了回调方法，用户自己改写借口即可