При составлении алгоритмов символического исчисления мы будем широко использовать объекты таких ранее созданных классов, как полиномы и векторы. Большинство данных и обслуживающих методов сосредоточим в заголовочном файле laplas.h (лучше дать в нем только интерфейсную часть, а реализационную перенести в файл laplas.cpp и включить его в состав программного проекта – последний вариант даст возможность использовать однократную предварительную компиляцию заголовочного файла при отладке программы, если не использовать в этом файле процедуры инициализации данных). Если в дальнейшем предполагается создать из файла реализации библиотечный файл, то разделение на интерфейсную и реализационную части обязательно.
Как всегда, в начало программы собираем все понадобившиеся при ее составлении заголовочные файлы интерфейса с библиотеками:
//Файл difequat.h
#include "vector.h"
#include "equation.h"
#include <windows.h>
#include <alloc.h>
#include "wlaplasm.rh"
#include "matrix.h"
/*Для упрощения записи конкретизируем типы векторных и матричных элементов */
typedef vector<double> dvector;
typedef matrix<double> dmatrix;
/*Шаблон структуры для сведений о корнях характеристического полинома*/
struct CA{
complex ROOT; //значение корня
complex A; //коэффициент для корня в оригинале
int J; //кратность
int O;}; //кому кратен
dmatrix SolMtr; /*Матрица решения и его производных*/
dvector userfunc; /*Вектор для произвольного возмущения*/
dvector tc; //Вектор аргумента решений
dmatrix FreqChMtr; /*Матрица для частотных характеристик*/
double w; //Частота среза
double Cufz; //Усиление на частоте среза
/*Класс обыкновенных линейных дифференциальных уравнений (ОЛДУ), содержащий методы работы, базирующиеся на операционном исчислении. */
class DifferentialEquation
{
//приватные данные
int rngl, rngr; /*Порядок уравнения - левая и правая часть*/
dvector nv; //Вектор начальных условий
double tend, tstep; /*Конечное время решения и шаг дискретизации по времени*/
int PointSolCnt; //Размерность вектора решения
//Полиномы левой, правой части и начальных условий
cpolynom PL,PR,PN;
/*Знаменатель и числитель изображения решения в виде полиномов*/
cpolynom PSOL,QSOL;
CA *R; /*указатель на структуру со сведениями о корнях*/
long RootCount; /*Количество корней общее - при наличии правой части уравнения оно не совпадает с его порядком слева*/
long rcount; /*Количество различных корней - за вычетом кратных*/
int cod; double omega, alpha; /*Сведения о стандартных возмущениях*/
public: //общедоступные члены класса
/*Конструктор должен получить коэффициенты уравнения, начальные условия, код стандартного возмущения и его параметры (частоту гармонической функции, постоянную экспоненты), конечное время решения и шаг*/
/*Инициализация полиномов с преобразованием коэффициентов в комплексные и реверсированием их последовательности - мы предполагаем, что вводятся коэффициенты и начальные условия начиная со старшей производной*/
int i;
for(i=0;i<(rngl+1);i++) PL[i]=complex(CFNL[i],0);
PL=PL.reverse();
for(i=0;i<(rngr+1);i++) PR[i]=complex(CFNR[i],0);
PR=PR.reverse();
//Создадим также вектор начальных условий
nv=dvector(rngl);
for(i=0;i<rngl;i++)
nv[i]=NU[rngl-i-1];
/*Формирование полинома начальных условий (числителя изображения решения при свободном движении. Если будет ненулевая правая часть, то ее изображение надо прибавить к результату вычисления PN)*/
//p-полином первой степени р+0
cpolynom p(2);
p[0]=complex(0.0,0.0);p[1]=complex(1.0,0.0);
cpolynom D=PL,tpn;
for( i=0;i<rngl;i++)
{
D=D/p;//делим полином на p, понижая степень на 1
/*прибавляем произведение значения производной i-го порядка в начальной точке на полиномиальный коэффициент*/
for(int j=0;j<rngl;j++) D[i]*=nv[i];
tpn=nv[i]*D;
PN+=tpn;
}
/*Сформируем числитель QSOL и знаменатель PSOL изображения решения - это полиномиальная дробь*/
cpolynom rone;
rone[0]=complex(1.0,0.0);/* веществ. полиномиальная единица */
//Если возмущения нет
if(cod==0) {PSOL=PL; QSOL=PN;}
/*Если это единичный импульс или произвольная функция*/
if(cod==1||cod==6) {PSOL=PL; QSOL=PR+PN;}
/*Если на входе ступенька - изображение 1/p
if(cod==2) {QSOL=PN*p+PR; PSOL=PL*p;}
/*Если на входе синусоида - изображение w/(p2+w2).*/
/*Если объект конструируется для расчета частотных характеристик*/
if(cod==7)
{
FreqChar(); /*Вызываем метод расчета частотных характеристик*/
RootCount=0;//Нет необходимости вычислять корни
return; //Досрочно покидаем конструктор
}
/*Определим теперь общее количество корней характеристического полинома*/
RootCount=PSOL.getm()-1;
//После определения количества корней выделим память для хранения сведений о них*/
R=(CA*)farcalloc(RootCount,sizeof(CA));
memset(R,0,RootCount*sizeof(CA));
//Конструируем хранители решений
double j;
tc=dvector(PointSolCnt); //Вектор аргумента
for(i=0,j=0.0;i<PointSolCnt;i++,j+=tstep)
tc[i]=j;
/*Матрица решений на RootCount строк - для функции и ее производных*/
SolMtr=dmatrix(RootCount+1,PointSolCnt);
/*Сразу занесем в матрицу решений начальные значения функции и ее производных*/
if(RootCount>0)
for(i=0;i<RootCount;i++) SolMtr[i][0]=nv[i];
GetRoot(); //Вычисляем корни х-го полинома
}//Конец конструктора
/*Функция для вычисления корней характеристического полинома*/
void DifferentialEquation::GetRoot()
{
int i,j;
//ищем корни характеристического полинома
cvector polyroot=newton(PSOL);
//Ограничим точность вычисления корней вблизи нуля
for(i=0;i<RootCount;i++)
{
if(fabs(real(polyroot[i]))<1e-7)
polyroot[i]=complex(0,imag(polyroot[i]));
if(fabs(imag(polyroot[i]))<1e-7)
polyroot[i]=complex(real(polyroot[i]),0);
}
//Заносим в структуру значения корней
for(i=0;i<RootCount;i++) R[i].ROOT=polyroot[i];
/*определяем кратность каждого корня и заносим в структуру; признаком -1 пометим корни, уже проверенные на кратность. J=1 будет корень первой кратности корню О и т.д.*/
int repeat;
rcount=RootCount;/* Вначале предполагаем, что все корни различны*/
for(i=0;i<RootCount;i++)
{
repeat=1;
if(R[i].J!=-1) //корень еще не встречался
{
R[i].J=1; //количество повторений i-го корня
for(j=i+1;j<RootCount;j++)
if((R[j].J!=-1)&&(R[i].ROOT==R[j].ROOT))
{
repeat++;
rcount--;R[j].J=-1;R[j].O=i;
/*устанавливаем признак того, что этот корень уже учтён*/
}
R[i].J=repeat;//кратность корня
}
}
}
/*Подпрограмма определения числителя изображения производной решения*/
cpolynom DifferentialEquation::Qsol(int der)
{
cpolynom Q=QSOL, p(2),sum;
//доформируем числитель изображения
p[0]=complex(0,0);p[1]=complex(1.0,0.0);
if(der>0 && cod!=6)
{
for(int i=0;i<der;i++)
sum+=PN[RootCount-i-1]*p^(der-i-1);
Q=Q*(p^der)-PSOL*sum;
}
return Q;
}
/*Функция, вычисляющая коэффициенты оригинала для всех корней; результат помещается в массив структур R*/
void DifferentialEquation::GetCoeffOrigin()
{
complex ap, tp;/*Для значений числителя и знаменателя*/
/*при подстановке значения корня числитель может оказаться полиномом или просто числом*/
int qi=QSOL.getm();//Выясняем это
int pi=PSOL.getm();
//Если это число
if(qi==1) ap=QSOL[0];/*Его и используем в качестве значения числителя*/
if(pi==1) tp=PSOL[0];
if(RootCount==1)//Единственный некратный корень
{
/*Вычисляем значение производной знаменателя изображения решения при подстановке в нее значения единственного корня*/
if(pi>1) tp=derive(PSOL,1)(R[0].ROOT);
if(qi>1) ap=QSOL(R[0].ROOT);/*Если числитель - полином, подставляем в него значение корня*/
R[0].A=ap/tp;
return;
}
//Если корней больше одного
long i,j,k,l;
cpolynom CH=QSOL,ZN=PSOL,RCH,RZN;
cpolynom pw=2;
cpolynom A;
//dcomplex ch;
for(k=0;k<RootCount;k++)
{
if(R[k].J==1) //Для некратного корня
{ //Формируем полином pw вида p+0
cpolynom pw(2);
pw[0]=-R[k].ROOT;pw[1]=complex(1.0,0.0);
if(qi>1) ap=QSOL(R[k].ROOT);
if(pi>1) tp=(PSOL/pw)(R[k].ROOT);
R[k].A=ap/tp;
}
//Если корень кратный
if(R[k].J>1)
{
cpolynom pw(2);pw[0]=-R[k].ROOT;
pw[1]=complex(1.0,0.0);
cpolynom pw1=(pw^R[k].J);
//До кратности этого корня
for(j=1;j<=R[k].J;)
{
if(j==1)
{
if(qi>1) ap=QSOL(R[k].ROOT);
if(pi>1) tp=(PSOL/pw1)(R[k].ROOT);
R[k].A=ap/tp;j++;
}
else
{
for(i=k+1;i<RootCount;i++)
{
//Если нашли кратный текущему
if(R[i].O==k)
{
cpolynom tmp=PSOL/pw1;
//Восстанавливаем знаменатель
cpolynom chisl=QSOL, znamen=PSOL/pw1,
dchisl, dznamen, m1, m2;
/*Дифференцируем дробь изображения
решения j-1 раз*/
for(l=1;l<=(j-1);l++)
{
dchisl=derive(chisl,1);
dznamen=derive(znamen,1);
m1=dchisl*znamen;
m2=dznamen*chisl;
chisl=m1-m2;
znamen^=2;
}
R[i].A=chisl(R[k].ROOT)/
(znamen(R[k].ROOT)*fact(j-1));
j++;
}
}
}
}
}
}
}
/*Подпрограмма, возвращающая значение функции решения ОЛДУ или ее производной заданного порядка при заданном t. Она реализует формулу
=f(t),
где
.*/
double DifferentialEquation::GetValue(double t)
{
int k,i,index;
double d;
complex result(0,0);
if(RootCount>0)
{
for(k=0;k<rcount;k++)
{
index=0;
if(R[k].J==1)//Если корень простой
result+=R[k].A*exp(R[k].ROOT*t);
if(R[k].J>1)//Если корень кратный
{
result+=R[k].A*exp(R[k].ROOT*t);index++;
for(i=k+1;i<rcount;i++)
{
if((R[i].O==k)&&(R[i].J==-1))
{
double dpw=R[k].J-index;
result+=R[k].A*pow(t,dpw)*
exp(R[k].ROOT*t)/fact(R[k].J-index);
index++;
}
}
}
}
d=real(result);
return d;
}
return 0;
}
/*Подпрограмма решения дифуравнения - она просто вызывает имеющиеся методы*/
void DifferentialEquation::Sol()
{
int i,k;
cpolynom Q=QSOL;
if(cod==6)
{
GetCoeffOrigin();
/*Вычисляем реакцию на импульс и пишем в 0-ю строку SolMtr*/
for(k=1;k<PointSolCnt;k++)
SolMtr[0][k]=GetValue(tc[k]);
/*Реакцию на произвольное возмущение разместим в 1-й строке SolMtr*/
for(i=1;i<PointSolCnt;i++)
{
SolMtr[1][i]=0;
for(k=0;k<i;k++)
SolMtr[1][i]+=SolMtr[0][i-k]*
userfunc[i]*tstep;
}
}
else
for(i=0;i<RootCount;i++)
{
QSOL=Qsol(i);
GetCoeffOrigin();
for(k=1;k<PointSolCnt;k++)
SolMtr[i][k]=GetValue(tc[k]);
QSOL=Q;
}
}
/* Подпрограмма вычисления частотных характеристик. Вначале подбирает частоту среза - то есть определяет верхнюю границу диапазона частот для исследования - по заданному пользователем коэффициенту снижения коэффициента усиления по отношению к нулевой частоте. Для этого частота меняется по простому алгоритму с переменным по величине и знаку шагом и вычисляется выходная амплитуда на этой частоте как модуль передаточной функции после подстановки в нее p=jw.*/
void DifferentialEquation::FreqChar()
{
FreqChMtr=dmatrix(5,PointSolCnt);
complex jone=complex(0.0,1.0);// мнимая единица
//QSOL=PR;PSOL=PL;
double KU0=fabs(real(PR[0]/PL[0])); /*Усиление на нулевой частоте*/
double KUZ=Cufz*KU0; //Заданное усиление
double wstep=1.0, //Начальный шаг по частоте
QW, //Значение числителя АФХ при частоте w.
PW, //Значение знаменателя АФХ при частоте w.
KU=KU0, /*Усиление при частоте w - вначале равно KU0.*/
RQ, //Вещественная часть значения числителя
IQ, //Мнимая часть значения числителя
RP, //Вещественная часть значения знаменателя
IP; //Мнимая часть значения знаменателя
if(PR.getm()==1)
QW=real(PR[0]); //Если справа - просто число
//Подбираем частоту среза
for(w=wstep;;w+=wstep)
{
if(PR.getm()>1) /*Если изображение правой части – полином*/
QW=sqrt(real(PR(w*jone))*real(PR(w*jone))+
imag(PR(w*jone))*imag(PR(w*jone)));
PW=sqrt(real(PL(w*jone))*real(PL(w*jone))+
imag(PL(w*jone))*imag(PL(w*jone)));
KU=QW/PW; //Усиление на текущей частоте
if((KU< 1.1*KUZ) && (KU>0.9*KUZ))
break; //Если попали в диапазон
if((KU<0.9*KUZ) && (wstep>0)) wstep=-0.1*wstep;
if((KU>1.1*KUZ) && (wstep<0)) wstep=-0.1*wstep;
}
/*Заполним вектор частоты в матрице - мы разместим его в 0-й строке*/
wstep=w/PointSolCnt;
int i;
for(i=0;i<PointSolCnt;i++)
FreqChMtr[0][i]=i*wstep;
//Вычисляем частотные характеристики
if(PR.getm()==1) {RQ=real(PR[0]); IQ=0;}
for(i=0;i<PointSolCnt;i++)
{
if(PR.getm()>1)
{
RQ=real(PR(jone*FreqChMtr[0][i]));
IQ=imag(PR(jone*FreqChMtr[0][i]));
}
RP=real(PL(jone*FreqChMtr[0][i]));
IP=imag(PL(jone*FreqChMtr[0][i]));
FreqChMtr[1][i]=(RQ*RP+IQ*IP)/(RP*RP+IP*IP);
//Вещественная частотная
FreqChMtr[2][i]=(IQ*RP-RQ*IP)/(RP*RP+IP*IP);
//Мнимая частотная
FreqChMtr[3][i]=sqrt(FreqChMtr[1][i]*
FreqChMtr[1][i]+FreqChMtr[2][i]*
FreqChMtr[2][i]); //Амплитудная
/*При расчете фазовой х-ки учтем скачек тангенса при фазе p/2*/
if(FreqChMtr[1][i]>0)
FreqChMtr[4][i]=atan(FreqChMtr[2][i]/
FreqChMtr[1][i]);
if(FreqChMtr[1][i]<=0)
FreqChMtr[4][i]=-M_PI+atan(FreqChMtr[2][i]/
FreqChMtr[1][i]);
//Если скачек не учитывать, то можно проще:
/*FreqChMtr[4][i]=atan(FreqChMtr[2][i]/
FreqChMtr[1][i]);//Фазовая */
}
}
/*Программа для решения ОЛДУ с постоянными коэффициентами символическим методом - файл difequat.cpp*/
#include "difequatm.h"
int sc; //Количество точек решения
char Userfile[80]; /*Для имени файла с произвольным возмущением*/
/*Для графического отображения функций решения объявим массив структур типа POINT, в который потом перенесем пересчитанные в координаты значения отображаемой функции и значения ее аргумента */
POINT *pt;
int DeriveNumber=0; //Порядок производной
/*Переменные для хранения размеров клиентской области окна вывода */
int cxClient, cyClient;
int Cod, //Код возмущения
Rngl,Rngr; //Порядок уравнения слева и справа
double Tend,Tstep, /*Конечное время решения и шаг по времени*/
Omega,Alpha; //Параметры гармоник и экспоненты
char* buf; //Буфер приема строк ввода пользователя
char*tmp; //Для приема указателя на слово от strtok
HINSTANCE hInst; /*Для сохранения идентификатора приложения*/
//Главная функция программы
#pragma argsused
int PASCAL WinMain(HINSTANCE hInstance, HINSTANCE hPrevInstance, LPSTR lpCmdLine, int nCmdShow)
{
hInst=hInstance;
WNDCLASSEX wc;
wc.cbSize=sizeof(wc);
wc.style=CS_HREDRAW|CS_VREDRAW;
wc.lpfnWndProc=(WNDPROC)WndProc;
wc.cbClsExtra=0;
wc.cbWndExtra=0;
wc.hInstance=hInstance;
wc.hIcon=(HICON)LoadIcon(0,IDI_WINLOGO);
wc.hCursor=(HCURSOR)LoadCursor(0,IDC_ARROW);
wc.hbrBackground=(HBRUSH)GetStockObject
(COLOR_BACKGROUND);//COLOR_BACKGROUND;
wc.lpszMenuName="Laplas";
wc.lpszClassName="Numerical Methods Class";
wc.hIconSm=0;
if(!RegisterClassEx(&wc))
{
MessageBox(0,"Не удалось зарегистрировать класс окна",0, MB_OK|MB_ICONEXCLAMATION);
return 0;
}
HWND hwnd=CreateWindowEx(WS_EX_CONTEXTHELP, "Numerical Methods Class", "Численные методы. Символический метод решения ОЛДУ",
WS_OVERLAPPEDWINDOW,
CW_USEDEFAULT,
CW_USEDEFAULT,
CW_USEDEFAULT,
CW_USEDEFAULT,
0,
0,
hInstance,
);
ShowWindow(hwnd,nCmdShow);
UpdateWindow(hwnd);
//userfile= new char[250];
MSG msg;
while(GetMessage(&msg,0,0,0))
{
TranslateMessage(&msg);
DispatchMessage(&msg);
}
delete[] pt;
return 0;
}
/* Массив значений отображаемой функции и ее аргумента нужно преобразовать в массив значений координат в окне отображения. Для такого пересчета напишем подпрограмму GetCoord - в качестве аргументов ей понадобятся массивы значений функции и аргумента. */
//Структура для смещений координатных осей
struct OFF{long x,y;}off;
void getCoord(dvector y, dvector x)
{
long i;
/*Определим наибольшее, наименьшее значение функции, диапазон изменения- хотя было бы неплохо получать эти значения из класса vector */
/* Можно нарисовать координатные оси - но мы этого не делаем, ограничившись качественной картиной процессов; предоставляем это вам для самостоятельной проработки. Один из вариантов:
=f(t),
.*/