global N Kn xL xR M0 imax dy alpha beta gamma Pr omega %clear all Eqmethod = 'NCCR'; M0=2; Kn = 1; L = 10*Kn; xL = -L; xR = L; N = 2; omega = 1; Pr = 0.661; % alpha = -0.31966;%-2/5; % beta = -0.281638;%-0.28; % gamma = 0.0638958;%0.5; alpha =-2/5; beta = 0;%-0.28; gamma = 0;%0.5; imax=1500; dy=(xR-xL)/(imax+1); Y = (xL+dy:dy:xR-dy); X0 = zeros(N*imax,1); for i = (1:imax) a1 = (sqrt(5/3)*(3 + 5*M0^2))/(8*M0); b1 = -(sqrt(15)*(-1 + M0^2)*coth(L))/(8*M0); X0((i-1)*N+1) = a1+b1*tanh(Y(i)); a2 = (-3 + 5*M0^2*(6 + M0^2))/(32*M0^2); b2 = ((-3 - 2*M0^2 + 5*M0^4)*coth(L))/(32*M0^2); X0((i-1)*N+2) = a2+b2*tanh(Y(i)); end % for i = (1:imax) % a1 = 0.4841229182759/M0 + 0.8068715304599*M0; % b1 = (0.4841229182759 - 0.4841229182759*M0^2)/(L*M0); % X0((i-1)*N+1) = a1+b1*(Y(i)); % % a2 = 0.9375 - 0.09375/M0^2 + 0.15625*M0^2; % b2 = (-0.09375 - 0.0625*M0^2 + 0.15625*M0^4)/(L*M0^2); % X0((i-1)*N+2) = a2+b2*(Y(i)); % % end %options = optimset('Display','iter','Jacobian','off', 'TolFun', 1.00e-010,'MaxFunEvals',200000, 'MaxIter', 10); options = optimoptions('fsolve','Display','iter','Jacobian','off', 'TolFun', 1e-10, 'MaxIter', 15); if(strcmpi(Eqmethod, 'nsf')) [X,exitflag] = fsolve(@NSF1DShockC,X0,options); %[X,exitflag] = newton(@NSF1DShockC,X0); %[X,exitflag] = Broyden(@NSF1DShockC,X0); elseif (strcmpi(Eqmethod, 'NCCR')) [X,exitflag] =fsolve(@CCR1DShockC,X0,options); elseif (strcmpi(Eqmethod, 'GRAD')) [X,exitflag] =fsolve(@Grad13Manuel,X0,options); end V = zeros(size(Y)); TEMP = zeros(size(Y)); for i = (1:imax) V(i ) = X((i-1)*N+1); TEMP(i ) = X((i-1)*N+2); end DEN = (sqrt(5./3).*M0)./V; SYY = 1 + (5.*M0.^2)./3 - (sqrt(5./3).*M0.*(V.^2 + TEMP))./V; QY = (5.*sqrt(15).*M0.^3 - 18.*V - 30.*M0.^2.*V + 3.*sqrt(15).*M0.*(5 + V.^2 - 3.*TEMP))./18; datafile=strcat(Eqmethod,'x1DShock','Kn', num2str(Kn),'Size',num2str(imax),'Ma',num2str(M0),'alpha',num2str(alpha),'beta',num2str(beta),'gamma',num2str(gamma),'.dat'); % name of the data file(s) where data is printed dlmwrite(datafile, [Y/Kn; (DEN);(TEMP); QY; V; TEMP], 'delimiter', '\t','precision', 10) fprintf('Name of DataFile is : %s\n',datafile); plot(Y/Kn,V,'r') % hold on % plot(Y/Kn,TEMP,'b') % plot(Y/Kn,DEN,'-m') % obervations: M0=2; % Kn = sqrt(pi/8)*5/4; % L = 10*Kn; %M0=2; %imax = 500: tiny wiggles before shock % 18 15019 1.14772e-025 1.4782e-007 8.47e-013 0.894 % decreasing L to 5*Kn;seems takes lot of time to conerse...stoped % Increasing L to 15*Kn;imax = 750% wigles are still here % increasinf imax = 1500 % eems they damped % increasinf imax = 2000 and L = 10*Kn % they damped % strating SNC from grad solution