%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Program to perform 4D-Var data assimilation on the % Lorenz equations % % % (c) 2002 Data Assimilation Research Centre % % Original Fortran program by Marek Wlasak % Converted to Matlab by Amos Lawless % % Change history: % Changes to make code more robust and add further options % (Amos Lawless) % 16/03/04 Corrections to code % Remove high memory allocations % Ensure windows appear in correct place % (Amos Lawless) % % List of main variables % a: sigma coefficient in equations % r: rho coefficient in equations % b: beta coefficient in equations % % ta: Length of assimilation window % tf: Length of forecast window % h: Time step for numerical scheme % freq: Frequency of observations in time steps % tstep: Number of time steps for assimilation % fcstep: Total number of time steps for assimilation + forecast % % [xx,yy,zz]: Truth values of trajectory % [xf,yf,zf]: First guess values of trajectory % [x,y,z]: Calculated values of trajectory and final analysis % [lx,ly,lz]: Gradient values calculated from adjoint % % [datx,daty,datz]: Observation values % D: Observation weighting matrix % max_iterations: Maximum number of minimization iterations % tolerance: Convergence tolerance for minimization % s: Calculated step length in minimization % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %--------------------------------------------------------------------- % 1. Set up %--------------------------------------------------------------------- % Close any open figure windows g=gcf; for i=1:g close(i) end clear % Determine screen size and set up figure positions ss=get(0,'ScreenSize'); fig_width=0.55*ss(3); fig_height=0.4*ss(4); pos_1=[0.02*ss(3),0.55*ss(4),fig_width,fig_height]; pos_2=[0.02*ss(3),0.05*ss(4),fig_width,fig_height]; pos_3=[0.6*ss(3),0.45*ss(4),0.4*ss(3),0.45*ss(4)]; pl=0.6*ss(3); pb=0.05*ss(4); pw=0.3*ss(3); ph=0.3*ss(4); pos_4=[pl,pb,pw,ph]; pos_5=[0.02*pw,0.8*ph,0.95*pw,0.2*ph]; pos_6=[0.02*pw,0.02*ph,0.95*pw,0.75*ph]; % format long %parameters a=1.0d1; % sigma b=8.0d0/3.0d0; % rho r=2.8d1; % beta % %--------------------------------------------------------------------- % 2. Input truth parameters and calculate truth %--------------------------------------------------------------------- % Inputs 1 xstr='True value of x at t=0'; ystr='True value of y at t=0'; zstr='True value of z at t=0'; true_str=inputdlg({xstr,ystr,zstr},'Initial data'); truex=str2num(true_str{1}); truey=str2num(true_str{2}); truez=str2num(true_str{3}); % Inputs 2 xstr='Length of assimilation window'; x2str='Length of subsequent forecast'; ystr='Time step'; zstr='Frequency of observations (in time steps)'; temp_str=inputdlg({xstr,x2str,ystr,zstr},'Run information'); ta=str2num(temp_str{1}); tf=str2num(temp_str{2}); h=str2num(temp_str{3}); freq=str2num(temp_str{4}); tstep = ta/h+1; fcstep = (ta+tf)/h; % Run truth guessx(1)=truex; guessy(1)=truey; guessz(1)=truez; j=0; [xx,yy,zz]=modeuler(fcstep,h,j,guessz,guessy,guessx,a,r,b); %--------------------------------------------------------------------- % 3. Calculate observations and plot truth and observations %--------------------------------------------------------------------- D=zeros(tstep,1); n_obs = (tstep-1) / freq + 1; datx=zeros(n_obs,1); daty=zeros(n_obs,1); datz=zeros(n_obs,1); l_noise=menu_asl('Noise on observations?','No','Yes'); l_noise=l_noise-1; if l_noise==1 sc_x_noise=randn(n_obs,1); sc_y_noise=randn(n_obs,1); sc_z_noise=randn(n_obs,1); sd_str=inputdlg('Variance of observation error'); sd=str2num(sd_str{1}); var=sqrt(sd); RX=var*sc_x_noise; RY=var*sc_y_noise; RZ=var*sc_z_noise; else for i=1:n_obs; RX(i)=0.0; RY(i)=0.0; RZ(i)=0.0; end end % Set up data and matrix D j=1; for i=1:freq:tstep datx(i)=xx(i)+RX(j); daty(i)=yy(i)+RY(j); datz(i)=zz(i)+RZ(j); D(i) = 1.; j=j+1; end % Plot truth and observations xvals=0:fcstep; vec=1:freq:tstep; v=vec-1; for i=vec x_ob(i)=datx(i); z_ob(i)=datz(i); end % x_ob and z_ob are temporary arrays for plotting % h1=figure('Position',pos_1); clf; subplot(2,1,1) plot(xvals,xx) hold on plot(v,x_ob(vec),'om') hold on xlabel('Time step') ylabel('x') title('Solution for x') legend('Truth','Observations') % subplot(2,1,2) plot(xvals,zz) hold on plot(v,z_ob(vec),'om') hold on xlabel('Time step') ylabel('z') title('Solution for z') legend('Truth','Observations') %--------------------------------------------------------------------- % 4. Input initial guess and run model from this %--------------------------------------------------------------------- xstr='Initial guess of x at t=0'; ystr='Initial guess of y at t=0'; zstr='Initial guess of z at t=0'; guess_str=inputdlg({xstr,ystr,zstr},'Initial guess'); guessx(2)=str2num(guess_str{1}); guessy(2)=str2num(guess_str{2}); guessz(2)=str2num(guess_str{3}); % j=1; % Plot first guess [xf,yf,zf]=modeuler(fcstep,h,j,guessz,guessy,guessx,a,r,b); subplot(2,1,1) plot(xvals,xf,':') legend('Truth','Observations','First guess') hold on subplot(2,1,2) plot(xvals,zf,':') legend('Truth','Observations','First guess') hold on %--------------------------------------------------------------------- % 5. Input minimization info and perform minimization %--------------------------------------------------------------------- xstr='Maximum number of iterations'; ystr='Tolerance'; y2str='Outer loops'; zstr={'30','1d-5','2'}; min_str=inputdlg({xstr,ystr,y2str},'Minimization control',1,zstr); max_iterations=str2num(min_str{1}); tolerance=str2num(min_str{2}); outer_loops=str2num(min_str{3}); % k=1; costplot=zeros(1); normplot=zeros(1); for i_outer=1:outer_loops %%%% Start of outer loop %%%% disp(['Outer loop ' num2str(i_outer)]) [x,y,z]=modeuler(tstep,h,i_outer,guessz,guessy,guessx,a,r,b); x0_p(1:2) = 0.0d0; y0_p(1:2) = 0.0d0; z0_p(1:2) = 0.0d0; % Calculate innovations for i=1:freq:tstep x_d(i) = x(i) - datx(i); y_d(i) = y(i) - daty(i); z_d(i) = z(i) - datz(i); end j=2; [f,g] = calcfg_inc(tstep,h,x,y,z,x_d,y_d,z_d,... x0_p(j),y0_p(j),z0_p(j),a,r,b,D,freq); cost(j)=f; lx=g(1); ly=g(2); lz=g(3); lnorm(j)=sqrt(lx(1)*lx(1)+ly(1)*ly(1)+lz(1)*lz(1)); % minimisation via least squares l_converged=0; if lnorm(j)tolerance & l_converged==0 % 70 s(j)=0.5d0; % normalise gradient vector gx=lx(1)/lnorm(j); gy=ly(1)/lnorm(j); gz=lz(1)/lnorm(j); x0_p(j+1) = x0_p(j)-s(j)*gx; y0_p(j+1) = y0_p(j)-s(j)*gy; z0_p(j+1) = z0_p(j)-s(j)*gz; [f,g] = calcfg_inc(tstep,h,x,y,z,x_d,y_d,z_d,... x0_p(j+1),y0_p(j+1),z0_p(j+1),a,r,b,D,freq); cost(j+1)=f; lx=g(1); ly=g(2); lz=g(3); lnorm(j+1)=sqrt(lx(1)*lx(1)+ly(1)*ly(1)+lz(1)*lz(1)); if (j-1>max_iterations) % Max number of iterations reached l_converged=1; disp('Failed to converge in maximum number of iterations') j=j-1; else Iteration=j-1 while (lnorm(j+1) >= lnorm(j) & s(j) > tolerance ... & l_converged == 0) % 75 s(j)=s(j)*0.5d0; x0_p(j+1) = x0_p(j)-s(j)*gx; y0_p(j+1) = y0_p(j)-s(j)*gy; z0_p(j+1) = z0_p(j)-s(j)*gz; [f,g] = calcfg_inc(tstep,h,x,y,z,x_d,y_d,z_d,... x0_p(j+1),y0_p(j+1),z0_p(j+1),a,r,b,D,freq); cost(j+1)=f; lx=g(1); ly=g(2); lz=g(3); lnorm(j+1)=sqrt(lx(1)*lx(1)+ly(1)*ly(1)+lz(1)*lz(1)); % Convergence test on norm of gradient if (abs((lnorm(j+1)-lnorm(j)))