Преобразование Лапласа

[lenachk@]

Здравствуйте.
Может у кого есть алгоритмы для обратного преобразования Лапласа?
Само изображение Лапласа имеет вид (1/(s+1)?

 

[????]

Ссылка скрыта от гостей

Ссылка скрыта от гостей

Ссылка скрыта от гостей

 

[lenachk@]

Спасибо (знаю в карман не спрячешь) , но вопрос все еще нерешен..
1.если кто может — прокоментируйте строки алгоритма.
как проверить правильность результата . есть прямое и обратное преобразование.


обратное :

#include "stdafx.h"

#include <iostream>
#include <cmath>
#include <complex>

using namespace std;
void InvLapTrans(complex<double> (*vs)(complex<double> s),double tstart,
double tend,double tdelta)
{
double t,vt;
int count;
static const complex<double> z[5] = {
complex<double>(11.83009373916819,1.593753005885813),
complex<double>(11.22085377939519,4.792964167565670),
complex<double>(9.933383722175002,8.033106334266296),
complex<double>(7.781146264464616,11.36889164904993),
complex<double>(4.234522494797000,14.95704378128156)}; 
static const complex<double> Kp[5] = {
complex<double>(16286.62368050479,-139074.7115516051),
complex<double>(-28178.11171305163,74357.58237274176),
complex<double>(14629.74025233142,-19181.80818501836),
complex<double>(-2870.418161032078,1674.109484084304),
complex<double>(132.1659412474876,17.47674798877164)};

t = tstart;
count = 0;
while (t <= tend) {
vt = 0.0;
for (int i=0;i<5;i++) {
vt = vt — real(vs(z[i] / t) * Kp[i]);
}
vt /= t;
cout<<count<<" "<<vt<<endl;
count++;
t += tdelta;
}
}

complex<double> vs(complex<double> s)
{

return 1.0/(s+1.0);
}

int main(int argc, char* argv[])
{
double tstart = 0.01;
double tend = 3.2;
double tdelta = 0.05;
InvLapTrans(vs,tstart,tend,tdelta);
return 0;
}

прямое преобразование:

// LaplInvTrans.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include <iostream>
#include <cmath>
#include <complex>

using namespace std;
void InvLaplace(int nt,double tr,int n,complex<double> p[],
complex<double> r[],double a[])
{
double t,dt,f;
int i,j;
complex<double> pt,ex;

t = 0.0;
dt = tr / nt;
for (i = 0;i< nt;i++) {
f = 0.0;
for (j=0;j<n;j++) {
if (imag(p[j]) > -1.0e-8) {
if (fabs(imag(p[j])) < 1.0e-8) {
/* contribution from real pole */ 
f += (real(r[j]) * exp(real(p[j])*t));

}
else {
pt = t * p[j];
ex = exp(pt);
ex *= r[j];
f += (2.0 * real(ex));

}

}
}
a[i] = f;

cout<<t<<"\t"<<a[i]<<endl;
t += dt;
}
}

Помогите. разобраться
есть еще один код , результаты аналог обратному

 using namespace std;
#include <complex>
#include <iostream>
typedef complex<double> cplxd;
const double PI = 3.141592653589793238;


cplxd F(const cplxd p)
{
cplxd tt;
cplxd tp(1,0); 
cplxd tf=tp/(p*p+tp*p);
tt=tf;			
return tt;
}

double invLaplace (double t, int n, double A) 
{
const int nmax = 500;
const double X=0.5*A/t, H=PI/t;
if (n >nmax) n = nmax;
// Euler summation (van Wijngaardenв algorithm)
double sum, sgn, wk[nmax+1];
int k = 0;
sgn = -1.0;
sum = 0.5*(wk[0]=0.5*real(F(cplxd(X,0))));
for (int m=1; m<n; m++, sgn=-sgn) 
{
double nxt, cur;
nxt = sgn*real(F(cplxd(X,m*H)));
for (int j=0; j<=k; j++) 
{
cur = wk[j]; wk[j] = nxt;
nxt = 0.5*(cur+nxt);
}
sum += (fabs(nxt)>fabs(cur)) ?
nxt : 0.5*(wk[++k]=nxt);  
}
return exp(0.5*A)/t * sum;
}

int main()
{
double t;
cplxd* p;
for (t=0.01;t<3.2;t+=0.05){
cout << "t= " << t << " Invlap=" << invLaplace(t,500,20) << endl;
}
system("PAUSE");
return 0;

}

а тут проблемы с математикой( откуда X=0.5*A/t параметр А и.все параметры.. )

Ссылка скрыта от гостей

такого нет B)

B)
Спасибо
2.

Ссылка скрыта от гостей

- как использовать (

Ссылка скрыта от гостей

для C++//