format short e; p.mut=0; p.muf=0; p.mub=0; p.dt=0.01; p.n=101; p.dx=1/(p.n-1); p.A=0.1; p.damp=0; try % % b =[ 10 9 8 7 6 5 4 %b=[3.76... % 3.1934 2 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.1] % b=[0.5 0.4 0.3 0.2 0.2 0.1 0.09 8.1312e-02 8.0900e-02 ... % 0.07 0.06 0.05 0.04 0.03]; %b=[3.1934 0.081312 0.0809 0.010371]; %b=1; %b=0.04; %omega=[1.87510407666,4.694091133,7.85475743824,10.9955407349, 14.137168391,17.278759532088237,20.420352251041251]; % b=4*pi^2./omega.^4; b=1; clear ZetaVals Err1Vals DisVals KappaVals x0Vals ThrustVals FricPowerVals TensionVals ClampPowerVals SpeedVals EnergyChangeVals; for count=1:length(b) p.B =b(count); n = p.n; dx=p.dx; p.Grid=linspace(0,1,n); p.tol = 1e-8; p.max_steps = 100; p.max_inner_ct = 5e2; p.t_final =200; p.Plot = 20; p.write=500; % % % filename1 = strcat('./','B.',... % num2str(p.B), '_mut.', num2str(p.mut),'_muf.',... % num2str(p.muf), '_mub.',num2str(p.mub),'_A.',num2str(p.A),'.mat'); % filename1 % save(filename1); outer_ct = 0; ZetaVals=[]; Err1Vals = []; DisVals=[]; x_guess=[]; KappaVals=[]; p.D1s=sparse(2:n-1,1:n-2,-1/2/dx*ones(n-2,1),n,n)+... sparse(2:n-1,3:n,1/2/dx*ones(n-2,1),n,n)+sparse([1;1;1],1:3,[-3/2/dx;2/dx;-1/2/dx],n,n)... +sparse([n;n;n],n-2:n,[1/2/dx;-2/dx;3/2/dx],n,n); p.D2s=sparse(1:n,1:n,[2/dx/dx;-2/dx/dx*ones(n-2,1);2/dx/dx],n,n)+... sparse(2:n-1,1:n-2,1/dx/dx*ones(n-2,1),n,n)+sparse(2:n-1,3:n,1/dx/dx*ones(n-2,1),n,n)... +sparse([1,1,1],2:4,1/dx/dx*[-5,4,-1],n,n)+sparse([n,n,n],[n-3:n-1],1/dx/dx*[-1,4,-5],n,n); p.Kappa=0.01*(p.Grid'-0.5*p.Grid'.^2-0.5); p.DtKappa=-0.1*(p.Grid'-0.5*p.Grid'.^2-0.5); p.Dtx0=-15; p.x0=0; p.theta = cumtrapz(p.Grid, p.Kappa); p.Dttheta=cumtrapz(p.Grid,p.DtKappa); p.DtZeta=cumtrapz(p.Grid,1i*p.Dttheta.*exp(1i*p.theta))+p.Dtx0+1i*2*pi*p.A; p.Zeta = cumtrapz(p.Grid, exp(1i*p.theta))+p.x0; % data=load('fB.2_mut.1_muf.0.01_mub.0.01_A.0.3.mat'); % KappaVals=data.KappaVals; % ZetaVals=data.ZetaVals; % p.Kappa=KappaVals(:,end); % p.DtKappa=(3/2*KappaVals(:,end)-2*KappaVals(:,end-1)+0.5*KappaVals(:,end-2))/p.dt; % p.theta = cumtrapz(p.Grid, p.Kappa); % % p.DtZeta=(3/2*ZetaVals(:,end)-2*ZetaVals(:,end-1)+0.5*ZetaVals(:,end-2))/p.dt; % % p.Zeta = ZetaVals(:,end); % p.Zeta = p.Zeta-real(p.Zeta(1,1)); % p.x0=real(p.Zeta(1,1)); % p.Dtx0=real(p.DtZeta(1,1)); % % % omega=[1.87510407666,4.694091133,7.85475743824,10.9955407349, 14.137168391,17.278759532088237,20.420352251041251]; % % p.y0=zeros(p.n,1); % p.Dty0=2*pi*p.A; % % for i=1:length(omega) % wi=omega(i); % lambdai=p.B*wi^4; % phi=(cosh(wi*p.Grid)-cos(wi*p.Grid)+(cosh(wi)+cos(wi))/(sin(wi)+sinh(wi))*(sin(wi*p.Grid)-sinh(wi*p.Grid)))'; % Norm=trapz(p.Grid,phi.*phi); % bi=4*pi^2*p.A*trapz(p.Grid,phi)/(-4*pi^2+lambdai)/Norm; % % p.Dty0=p.Dty0+2*pi*bi*phi; % end % p.Zeta=p.Grid'+1i*p.y0; % p.DtZeta=1i*p.Dty0; % % p.Kappa=(p.D2s*p.y0)./(1+(p.D1s*p.y0).^2).^1.5; % p.theta = cumtrapz(p.Grid, p.Kappa); p.t=p.dt; p.Theta0=0; p.Initial(2:p.n+1,1)=p.Kappa; p.Initial(1,1)=p.x0; [r]=Initial(p.Initial,p); KappaVals(:,1)=r.Kappa; ZetaVals(:,1)=r.Zeta; DisVals(1)=imag(r.Zeta(p.n,1)); x0Vals(1)=r.x0; ThrustVals(1)=r.thrust; FricPowerVals(1)=r.fricpower; % EnergyVals(1)=trapz(p.Grid,0.5*abs(r.DtZeta).^2+0.5*p.B*r.Kappa.^2); % ClampPowerVals(1)=p.B*(-3*r.Kappa(1)+4*r.Kappa(2)-r.Kappa(3))/2/p.dx*2*pi*p.A*cos(2*pi*p.t); TensionVals(:,1)=r.T; DtZetaVals(:,1)=r.DtZeta; DttZetaVals(:,1)=r.DttZeta; p.Zeta_2=r.Zeta; p.Zeta_1=p.Zeta; p.Zeta_0=r.Zeta1; p.x0_2=r.x0; p.x0_1=p.x0; p.x0_0=r.x01; x_guess=[]; ZetaStart=p.Zeta_1; KappaStart=p.Kappa; p.flag=0; k1=2; while (p.t < p.t_final) if p.flag==0 p.mut=10000; %p.muf=0.03*p.mut; %p.mub=0.03*p.mut; p.muf=0.01; p.mub=0.01; p.damp=0; p.flag=1; filename1 = strcat('./','fB.',... num2str(p.B), '_mut.', num2str(p.mut),'_muf.',... num2str(p.muf), '_mub.',num2str(p.mub),'_A.',num2str(p.A),'.mat'); filename1 %save(filename1); end p.t = p.dt*k1; if k1==1 x_guess(1,1)=r.x0; x_guess(2:p.n+1,1)=KappaStart; elseif k1==2 x_guess(2:p.n+1,1)=KappaVals(:,k1-1); x_guess(1,1)=x0Vals(k1-1); else x_guess(2:p.n+1,1) = 2*KappaVals(:,k1-1)-KappaVals(:,k1-2); x_guess(1,1)=2*x0Vals(k1-1)-x0Vals(k1-2); %%%%%%%%%%%% interpolation from previous 2 step solutions end %%%%%%%%%% %%Broyden Method using S-M Formula %%%%%%%%%%%%%%%%%%% if k1==2 A = GetJacobianInv(x_guess, p); end [r] = Fcn(x_guess,p); err = norm(r.out); err1 = err; inner_ct = 0; x_new = x_guess; while (err > p.tol) && (inner_ct < p.max_inner_ct) inner_ct = inner_ct + 1; smaller = 0; j = 0; f0 = r.out; d = -(A*r.out); for j = 0:p.max_steps-1 x_new = x_guess + (2^-j)*d; [r]=Fcn(x_new,p); err1 = norm(r.out); if (err1 < err) break; end end w = r.out - f0; q = x_new - x_guess; qA = (q'*A); denom = qA*w; A_update = (A*w - q)*qA/denom; %%%%%%%%%%Sherman-Morrison Formula%%%%%%%%%%% A = A - A_update; if mod(inner_ct, 20) == 0 [1 err1]; A = GetJacobianInv(x_new, p); inner_ct = inner_ct + 5; end if abs(denom) < 1e-15 || norm(q) < 1e-15 || ... norm(A_update) > 1e8 [2 denom norm(q) norm(A_update) err1]; A = GetJacobianInv(x_new, p); inner_ct = inner_ct + 20; end err = err1; x_guess = x_new; end outer_ct = outer_ct+1; ZetaVals(:,k1)=r.Zeta; DisVals(k1)=imag(r.Zeta(p.n,1)); KappaVals(:,k1)=r.Kappa; x0Vals(k1)=r.x0; ThrustVals(k1)=r.thrust; FricPowerVals(k1)=r.fricpower; TensionVals(:,k1)=r.T; DtZetaVals(:,k1)=r.DtZeta; DttZetaVals(:,k1)=r.DttZeta; % EnergyVals(k1)=trapz(p.Grid,0.5*abs(r.DtZeta).^2+0.5*p.B*r.Kappa.^2); SpeedVals(k1-1)=trapz(p.Grid,real(r.DtZeta)); if k1==2 dtKappa=r.Kappa-KappaVals(:,k1-1); % EnergyChangeVals(k1-1)=(EnergyVals(k1)-EnergyVals(k1-1))/p.dt; else dtKappa=(3/2*r.Kappa-2*KappaVals(:,k1-1)+0.5*KappaVals(:,k1-2))/p.dt; % EnergyChangeVals(k1-1)=(3/2*EnergyVals(k1)-2*EnergyVals(k1-1)+0.5*EnergyVals(k1-2))/p.dt; end EnergyChangeVals(k1-1)=trapz(p.Grid,real(r.DtZeta.*conj(r.DttZeta))+p.B*r.Kappa.*dtKappa); ClampPowerVals(k1-1)=p.B*(-3*r.Kappa(1)+4*r.Kappa(2)-r.Kappa(3))/2/p.dx*imag(r.DtZeta(1)); % if mod(k1,100)==0 % % r1=norm(ZetaStart-r.Zeta); % r2=norm(KappaStart-r.Kappa); % if(r1<1e-4&&r2<1e-4) % save(filename1,'ZetaVals', 'DisVals','KappaVals','ThrustVals'); % break; % else % % ZetaStart=r.Zeta; % KappaStart=r.Kappa; % end % end if mod(k1,p.Plot) == 0; time=p.dt:p.dt:p.t; figure(1); plot(real(r.Zeta), imag(r.Zeta), 'r-'); hold on; axis equal; figure(2); plot(time,DisVals,'m-'); end p.Zeta_0 = p.Zeta_1; p.Zeta_1 = p.Zeta_2; p.Zeta_2 = r.Zeta; p.x0_0 = p.x0_1; p.x0_1 = p.x0_2; p.x0_2 = r.x0; k1=k1+1; if p.t>5 && mod(k1,p.write)==0 save(filename1,'ZetaVals', 'DisVals','KappaVals','ThrustVals','FricPowerVals','ClampPowerVals','SpeedVals', 'EnergyChangeVals','TensionVals'); end end end catch ME ErrMsg = ME; save(filename1,'ZetaVals', 'DisVals', 'ErrMsg','KappaVals','ThrustVals','FricPowerVals','ClampPowerVals','SpeedVals', 'EnergyChangeVals'); end %save(filename1,'ZetaVals', 'DisVals','KappaVals');