Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Integration RUNGE_KUTTA in complex plane.

  1. Dec 14, 2009 #1
    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
     
  2. jcsd
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook

Can you offer guidance or do you also need help?
Draft saved Draft deleted