Решение краевой задачи для линейного дифференциального уравнения второго порядка методом конечных разностей на C++


При численном решении этой задачи методом конечных разностей отрезок [a,b] разбивают на равные части с шагом h точками xi, где h=(b−a)/n.

Заменяя производные правыми односторонними конечно — разностными отношениями для внутренних точек xi и концевых точек отрезка [a,b], и выполняя алгебраические преобразования, получим систему линейных алгебраических уравнений




Реализовать решение на C++ можно следующим образом:


//============================================================================
// Name        : mkr.cpp
// Author      : 
// Version     :
// Copyright   : Your copyright notice
// Description :
//============================================================================

#include <iostream>
#include <cmath>
using namespace std;

double p(double x){
	return exp(x);
}

double q(double x){
	return 2*x;
}

double f(double x){
	return pow(x,3);
}

void gauss(double** A, double* B, double* X, int n);

int main() {
	double alpha0=2; double alpha1=-2.5; double Ac=0;
	double beta0=3; double beta1=-3.4; double Bc=5;
	double a0=0.1; double b0=1.3;
	int n=4;
	//matrix A[n+1][n+1]
	double **A = new double*[n+1];
	for(int i=0; i<n+1; i++){
		A[i]=new double[n+1];
	}
	//===
	double *B = new double[n+1]; //vector B[n+1]
	double *X = new double[n+1]; //vector X[n+1]
	//Grid
	double h=(b0-a0)/n;
	for(int i=0; i<=n; i++){
		X[i]=a0+i*h;
		//cout << X[i] << endl;
	}
	//===
	cout << "h=" << h << endl;
	//calculate matrix A, B
	for(int i=0; i<=n-2; i++){
    	A[i][i]=h*h*q(X[i])-h*p(X[i])+1;
    	A[i][i+1]=h*p(X[i])-2;
    	A[i][i+2]=1;
    	B[i]=h*h*f(X[i]);
	}
	A[n-1][0]=alpha0*h-alpha1;
	A[n-1][1]=alpha1;
	A[n][n-1]=-beta1;
	A[n][n]=beta0*h+beta1;
	B[n-1]=h*Ac;
	B[n]=h*Bc;
	//===
	//print A
	for(int i=0; i<=n; i++){
		for(int j=0; j<=n; j++){
			cout << "A["<<i<<"]"<<"["<<j<<"]="<<A[i][j] <<" ";
		}
		cout << endl;
	}
	//===
	cout << endl;
	//print B
	for(int i=0; i<=n; i++){
		cout << "B["<<i<<"]="<<B[i] <<" ";
	}
	cout << endl;
	cout << endl;
	//===
	//solve A*X1=B
	double *X1 = new double[n+1]; //vector X1[n+1]
	gauss(A,B,X1,n+1); //solve
	//print X1
	for(int i=0; i<=n; i++){
		cout << "X1[" << i << "]=" << X1[i] << " ";
	}
	//===
	return 0;
}

void gauss(double** A, double* B, double* X, int n){
	int m=n+1;
	//create C[n][n+1]
	double **C = new double*[n];
	for(int i=0; i<n; i++){
		C[i]=new double[n+1];
	}
	//===
	//split A and B to C
	for(int i=0;i<n;i++){
		for(int j=0;j<n;j++){
			C[i][j]=A[i][j];
		}
		C[i][n]=B[i];
	}
	//===
	/*
	//print C
	for(int i=0;i<n;i++){
		for(int j=0;j<n+1;j++){
			cout << C[i][j] << " ";
		}
		cout << endl;
	}
	//===
	*/

	//forward
	for (int k=0; k<n-1; k++){
	    for (int i=k+1;i<n;i++){
	        for (int j=m-1; j>=k; j--){
		        C[i][j]=C[i][j]-C[i][k]*C[k][j]/C[k][k];
	        }
    	}
	}
	//reverse
	X[n-1]=A[n-1][m-2]/A[n-1][m-2];

	for (int i=n-2; i>=0; i--){
		double s=0;
		for (int j=i+1;j<m-1;j++){
			s=s+C[i][j]*X[j];
		}
		X[i]=(C[i][m-1] - s)/C[i][i];
	}

}
  • avatar
  • 0

0 комментариев

Только зарегистрированные и авторизованные пользователи могут оставлять комментарии.