About this blog

I feel this blog as a reflection of my thoughts to myself , and sometimes as a public diary, and the is my only friend to share my thoughts who says never a "oh no! ,you shouldn't....That is boring...."

fitting an exponentially decaying co-sine using differential equation

Octave or Matlab code for 
Sample function trial for fitting an exact data from exponentially decaying sine. 




##code starts from here
version:1 
x=-4:0.1:4;

u=exp(-2*x).*sin(pi*x+1);

du=exp(-2*x).*(pi*cos(pi*x+1)-2*sin(pi*x+1));
ddu=exp(-2*x).*(-4*pi*cos(pi*x+1)+(4-pi^2)*sin(pi*x+1));
beta=[-sum(du.*ddu); -sum(u.*ddu)];
A=[sum(du.*du), sum(u.*du);...
   sum(u.*du), sum(u.*u)];
param=A\beta;
b=param(1); c=param(2);
discrim=sqrt(b^2-4*c);
mu=(-b+discrim)/2;
nu=(-b-discrim)/2;
X=exp(mu.*x);
Y=exp(nu.*x);
beta=[sum(X.*u); sum(Y.*u)];
A=[sum(X.*X) sum(X.*Y);...
   sum(X.*Y) sum(Y.*Y)];
param=A\beta;
c=param(1); d=param(2);
v=real(c.*X+d.*Y);
plot(u-v); pause;  ##gives the pointwise error between original and fitted data.

Version:2 (with additive noise and central difference formula)
x=1:0.01:3;
dx=0.01;
noise=0; #1e-2*rand(size(x));
u=exp(-2*x).*sin(pi*x+1)+noise;
#u=200*exp(-2*x)+exp(3*x);
#du=exp(-2*x).*(pi*cos(pi*x+1)-2*sin(pi*x+1))+noise;
#ddu=exp(-2*x).*(-4*pi*cos(pi*x+1)+(4-pi^2)*sin(pi*x+1))+noise;

#du=diff(u); du=[2*du(1)-du(2) du];
du=u(3:end)-u(1:end-2); du=[2*du(1)-du(2) du 2*du(end)-du(end-1)];
du=du/dx/2;
ddu=du(3:end)-du(1:end-2); ddu=[2*ddu(1)-ddu(2) ddu 2*ddu(end)-ddu(end-1)];
#ddu=diff(du) ; ddu=[ddu 2*ddu(end)-ddu(end-1)];
ddu=ddu/dx/2;

######################################
##Begins Normal equations solving

beta=[-sum(du.*ddu); -sum(u.*ddu)];
A=[sum(du.*du), sum(u.*du);...
   sum(u.*du), sum(u.*u)];
param=A\beta;
b=param(1); c=param(2);
discrim=sqrt(b^2-4*c);
mu=(-b+discrim)/2;
nu=(-b-discrim)/2;
X=exp(mu.*x);
Y=exp(nu.*x);
beta=[sum(X.*u); sum(Y.*u)];
A=[sum(X.*X) sum(X.*Y);...
   sum(X.*Y) sum(Y.*Y)];
param=A\beta;
c=param(1); d=param(2);
v=real(c.*X+d.*Y);
#plot(x,v,'m');pause
semilogy(x,abs(u-v)./(abs(u)+eps),'b;Relative error;'); hold on;
semilogy(x,abs(u-v),'r;absolute error;');
title("blue:relatice red:absolute")
pause;



















కామెంట్‌లు లేవు: