Расчёт дифференциальных уравнений методом Эйлера в C++.

В статье "расчёт в маткаде переходных процессов в ёмкостном фильтре" рассматривался расчёт схемы изображенной на рисунке 1 в которой источник питания e(t) имитировал выходное напряжение однополупериодного однофазного выпрямителя. Расчёт проводился в программе MathCAD. Там же была составлена математическая модель схемы (рисунок 1) представляющая собой дифференциальное уравнение:

Решением этого уравнения (уравнения (1)) при заданных начальных условиях (начальном t и начальном Uc) является функция описывающая изменение напряжения конденсатора во времени.

Рисунок 2 - Схема фильтра с источником и обозначениями для расчёта

Ранее расчёт уравнения 1 был проведен в маткаде где уже присутствовала функция расчёта дифференциальных уравнений. Рассмотрим расчёт этого же уравнения реализованный на языке программирования C++. Исходный код программы приведен ниже:


#include <iostream>
#include <conio.h>
#include <math.h>
#include <fstream.h> 

using namespace std;

const double C=0.01;//ёмкость конденсатора C
const int Re=2;//внутреннее сопротивление источника Re
const int Rn=100;//сопротивление нагрузки Rn
//источник e(t)
double e(double t)
{
  double s;
  s=2*sin(1000*t);
  if (s>0)return s;
  else return 0;
}
//уравнение (1)
double sys(double t, double Uc)
{
  double dUc;
  dUc=(e(t)/(C*Re))-(Uc/(C*Re))-(Uc/(C*Rn));
  return dUc;    
}

int main ()
{
    double t[100];
    double Uc[100];
    double h=0.0005;//шаг интегрирования
   //начальные условия
    t[0]=0.0;
    Uc[0]=0.0;
    //открытие файлов для записи
    ofstream time("time.TXT");
    ofstream voltage("voltage.TXT");
    //Расчёт  методом Эйлера и вывод результатов
    for(int i=0;i<100;i++)
    {
       t[i]=t[0]+i*h;
       Uc[i+1]=Uc[i]+h*sys(t[i],Uc[i]);
       cout<<"t="<<t[i]<<"        Uc="<<Uc[i]<<endl;
       time<<t[i]<<endl;
       voltage<<Uc[i]<<endl;            
    }
    time.close();//закрытие
    voltage.close();//файлов
 
    _getch();
    return 0;
}

Параметры элементов схемы (1) заданны в виде констант: C, Re, Rn. Функция e принимает значение времени t и возвращает напряжение на входе однополупериодного однофазного выпрямителя. Функция sys принимает значения: времени t, напряжения конденсатора Uc и возвращает значение dUc/dt уравнения (1) при заданных Uc и t. В строчке:

 t[i]=t[0]+i*h; 

производится расчёт времени, где h - шаг интегрирования, в методе Эйлера выбирается достаточно малым для обеспечения требуемой точности. В строчке:

Uc[i+1]=Uc[i]+h*sys(t[i],Uc[i]);

рассчитывается значение напряжения конденсатора в следующий момент времени Uc[i+1] исходя из напряжения конденсатора в текущий момент времени Uc[i] и текущего момента времени t[i]. Рассчитанные значения, на каждой итерации, записываются в массивы, выводятся на экран и в файлы time.TXT и voltage.TXT.

  На рисунке 2, красной линией, показан график напряжения построенный в маткаде по значениям рассчитанным в программе и выведенным в файл. Чтобы создать в маткаде матрицу с выведенными в файл значениями нужно: 1) открыть файл (значения должны быть числовыми и располагаться каждый в своей строчке), 2) выделить значения, 3) скопировать (например нажав ctrl+c), 4) перейти в окно маткада и 5) вставить скопированные значения (например нажав ctrl+v) .

Рисунок 2 - График напряжения конденсатора Uc (красный) и источника e(t) (синий).

Из графика видно что форма напряжения, рассчитанная методом Эйлера, мало чем отличается (если отличается вообще) от полученной при расчёте в программе MathCAD. Также для вывода графика можно воспользоваться программой исходный код которой приведен в статье "программа для вывода графика на языке C++".