PDA

View Full Version : ntegration RUNGE_KUTTA in complex plane.

nit
12-14-2009, 04:48 AM
Hi,
I have one serious problem while solving rayleigh equation using blasius profile, in which so as to remove the singularity the intergration contour is defined in a complex plane.
4th order runge kutta is used but if the step size( h) is is in complex, it is giving some error. would anyone please help me solving it?
here is the matlab code. thank you...reply as soon as possible.

nstep=100;
pi=3.14;
error=1e-4;
ycr=1.35;
yci=1.35;
ymax=8;

nstep2=5*nstep-1;
ncirc=nstep;
ncirc2=3*ncirc+1;

for j=1:ncirc2
theta=pi*(j-1)/(ncirc2-1);
yr=ycr*(1-cos(theta));
yi=-yci*sin(theta);
Y(j)= (yr+i*yi);
end

for j=ncirc2:nstep2
yr=ycr*2+ymax*(j-ncirc2)/(nstep2-ncirc2);
Y(j)=(yr+i*0);
end

%y=zeros(1,61);
%z=zeros(1,61);
%p=zeros(1,61);
%u=zeros(1,61);

F(1)=0+0*i;
U(1)=0+0*i;
UP(1)=1.660287+0*i;
UPP(1)=0+0*i;
[y,z,p,u]=blasius1(nstep2,Y,F,U,UP,UPP);
pnew=p(1);
znew=z(nstep2);

UP(1)=1.05*UP(1);
[y,z,p,u]=blasius1(nstep2,Y,F,U,UP,UPP);
pold=p(1);
zold=z(nstep2);

residue=1e30;

while(residue>error)
U(1),UPP(1),F(1),Y,nstep2;
UP(1)=pnew-(((znew)-1)*((pold)-(pnew)))/((zold)-(znew));
[y,z,p,u]=blasius1(nstep2,Y,F,U,UP,UPP);
zold=znew;
pold=pnew;
znew=z(nstep2);
pnew=p(1);
residue=abs((z(nstep2))-1.0);
end
%plot(real(Y),real(z));
hold on
plot(real(Y),abs(p));
hold on
%plot(real(Y),real(y));
hold on
%plot(real(Y),real(u));
================================================== ===============
blasius.m(separate file):
function[y,z,p,u]=blasius1(num,Y,F,U,UP,UPP)

y(1)=F(1);
z(1)=U(1);
p(1)=UP(1);
u(1)=UPP(1);

for j=1num-1)

h=((Y(j+1))-(Y(j)))*0.5;
h=real(h);

A11=(h)*(z(j));
A12=(h)*(p(j));
A13=(-1*(h)*((y(j))*(p(j))));
A14=(-(1.4806)*(h)*(((y(j))*(u(j))+(z(j))*(p(j)))));

A21=(h)*((z(j))+(A12)*0.5);
A22=(h)*((p(j))+(A13)*0.5);
A23=(-1*(h)*(((y(j))+(A11)*0.5)*((p(j))+(A13)*0.5)));
A24=(-(1.4806)*(h)*((((y(j))+(A11)*0.5)*((u(j))+(A14)*0. 5)+((z(j))+(A12)*0.5)*((p(j))+(A13)*0.5))));

A31=(h)*((z(j))+(A22)*0.5);
A32=(h)*((p(j))+(A23)*0.5);
A33=(-1*(h)*(((y(j))+(A21)*0.5)*((p(j))+(A23)*0.5)));
A34=(-(1.4806)*(h)*((((y(j))+(A21)*0.5)*((u(j))+(A24)*0. 5)+((z(j))+(A22)*0.5)*((p(j))+(A23)*0.5))));

A41=(h)*((z(j))+(A32));
A42=(h)*((p(j))+(A33));
A43=(-1*(h)*(((y(j))+(A31))*((p(j))+(A33))));
A44=(-(1.4806)*(h)*((((y(j))+(A31))*((u(j))+(A34))+((z(j ))+(A32))*((p(j))+(A33)))));

y(j+1)=((y(j))+(((A11)+2*(A21)+2*(A31)+(A41))/6));
z(j+1)=((z(j))+(((A12)+2*(A22)+2*(A32)+(A42))/6));
p(j+1)=((p(j))+(((A13)+2*(A23)+2*(A33)+(A43))/6));
u(j+1)=((u(j))+(((A14)+2*(A24)+2*(A34)+(A44))/6));

end