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

Тема в разделе "C и С++ FAQ", создана пользователем lenachk@, 18 окт 2006.

Статус темы:
Закрыта.
  1. lenachk@

    lenachk@ Гость

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

    lenachk@ Гость

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


    обратное :
    Код (Text):
    #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;
    }

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


    Код (Text):
    // 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;
    }
    }
    Помогите. разобраться
    есть еще один код , результаты аналог обратному
    Код (Text):
     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. http://google.com/codesearch?q=Inverse+Lap...tnG=Search+Code - как использовать (
    http://www.google.com/codesearch?q=+Invers.../ilaplace.m#a0) для C++//
     
Загрузка...
Похожие Темы - Преобразование Лапласа
  1. WolfEater
    Ответов:
    0
    Просмотров:
    777
  2. jager
    Ответов:
    1
    Просмотров:
    1.295
  3. PahaStar
    Ответов:
    0
    Просмотров:
    985
  4. PahaStar
    Ответов:
    0
    Просмотров:
    845
  5. ask40
    Ответов:
    0
    Просмотров:
    966
Статус темы:
Закрыта.

Поделиться этой страницей