|
From: | Leo Baumann |
Subject: | Re: A problem with ilaplace() |
Date: | Tue, 14 Nov 2017 19:15:54 +0100 |
User-agent: | Mozilla/5.0 (Windows NT 10.0; WOW64; rv:52.0) Gecko/20100101 Thunderbird/52.4.0 |
Am 14.11.2017 um 18:22 schrieb Doug Stewart:
At the moment my script makes no use. I was preparing a general usable script to evaluate step function responses for electric networks. I could do that numeric with following Octave function, but the symbolic answer is much more interesting. - # Talbot suggested that the Bromwich line be deformed into a contour that begins # and ends in the left half plane, i.e., z ? -8 at both ends. # Due to the exponential factor the integrand decays rapidly # on such a contour. In such situations the trapezoidal rule converge # extraordinarily rapidly. # For example here we compute the inverse transform of F(s) = 1/(s+1) at t = 1 #----------------------------------------------------------------------------- #>> pkg load symbolic #>> syms s #>> F=1/(s+1) #F = (sym) # # 1 # ----- # s + 1 # #>> error=talbot(function_handle(F),1,24)-exp(-1) #ans = 1.6098e-015 #------------------------------------------------------------------------------- # # Talbot method is very powerful here we see an error of 1.61e-015 # with only 24 function evaluations # # Created by Fernando Damian Nieuwveldt # email:address@hidden # Date : 25 October 2009 # # Reference # L.N.Trefethen, J.A.C.Weideman, and T.Schmelzer. Talbot quadratures # and rational approximations. BIT. Numerical Mathematics, # 46(3):653 670, 2006. function ilt = talbot(f_s, t, N) h=2*pi/N; # Shift contour to the right in case there is a pole on the positive real axis : Note the contour will # not be optimal since it was originally devoloped for function with # singularities on the negative real axis # For example take F(s) = 1/(s-1), it has a pole at s = 1, the contour needs to be shifted with one # unit, i.e shift = 1. But in the test example no shifting is necessary shift=0; ans=0; for k=0:N theta=-pi+(k+1/2)*h; z=shift+N/t*(0.5017*theta*cot(0.6407*theta)-0.6122+0.2645i*theta); dz=N/t*(-0.5017*0.6407*theta/sin(0.6407*theta)^2+0.5017*cot(0.6407*theta)+0.2645i); ans=ans+exp(z*t)*f_s(z)*dz; endfor ilt=real((h/(2i*pi))*ans); endfunction Regards - Leo |
[Prev in Thread] | Current Thread | [Next in Thread] |