From e3a4af8e2c1e881e3ea8d39071ae38b42f1de958 Mon Sep 17 00:00:00 2001 From: Yang Liu Date: Tue, 25 Apr 2023 03:01:42 -0400 Subject: [PATCH] delete wrong files --- EMT/DT.m | 81 ---------- EMT/multiwindow.m | 375 ---------------------------------------------- 2 files changed, 456 deletions(-) delete mode 100644 EMT/DT.m delete mode 100644 EMT/multiwindow.m diff --git a/EMT/DT.m b/EMT/DT.m deleted file mode 100644 index 03e7b50..0000000 --- a/EMT/DT.m +++ /dev/null @@ -1,81 +0,0 @@ -function [t_new,x_new,y_new,LTE,LTE2,LTE3] = DT(t,x,f,GLO,h,K) -% DT performs one-step integration of differential transformation -% method to solve differential equation dx/dt = f(t,x,GLO); -% -% Input -% k current order, 0,1,2,... -% x current states -% f differential equation -% GLO parameters in the differential equation -% h time step -% K The degree of power series: X(0)+X(1)*h+...X(K)*h^K -% -% Output -% t_new next time instant -% x_new next states of state variables -% y_new next states of algebraic variables -% LTE local truncation error -% -% See also rk4, Parareal -% - - - %% Function handle to compute intermidate variables - fun_y = str2func(strcat(func2str(f),'_AE')); - fun_idx = str2func(strcat(func2str(f),'_idx')); - - %% Function handle for DT model X(k+1)=FB(X(0:k)) of differential equation dx/dt = f(x) - F = str2func(strcat(func2str(f),'_DT')); - - %% compute intermidate variables y - idx = fun_idx(GLO); - y = fun_y(x,idx,GLO); - - %% Define matrix X Y to store coefficients, Xh to store terms - Nx = length(x); - X = zeros(Nx,K+1); %%why NAN - - Ny = length(y); - Y = zeros(Ny,K+1); - - Xh= zeros(Nx,K+1); - Yh= zeros(Ny,K+1); - %% Initialize the 0th order coefficients - X(:,1) = x; - Y(:,1) = y; - - Xh(:,1) = x; % Term(0) = X(0)*h^0 - Yh(:,1) = y; - %% Perform recursive calculation of high order coefficients and terms - for k = 0:K-1 -% for k = 0:K-2 - %% DT Equation: X(k+1)=F(X(0:k)) or XX(k+2)=F(XX(1:k+1)) - [X,Y] = F(k,X,Y,GLO,idx); - Xh(:,k+2) = X(:,k+2)*h^(k+1); % Term(k+1) = X(k+1)*h^(k+1) or XX(k+2)*h^(k+1) - Yh(:,k+2) = Y(:,k+2)*h^(k+1); % Term(k+1) = X(k+1)*h^(k+1) or XX(k+2)*h^(k+1) - - end - - %% sum over all columns of matrix Xh to obtain x_new as a column vector - x_new = sum(Xh(:,1:end-1),2); - - %% -% y_new = sum(Yh,2); - y_new = fun_y(x_new,idx,GLO); - %% Update time - t_new = t + h; - eta = 0.0000000000000000001; - %% Local truncation error: last term; maximum value over all variables - % khuang: we can use relative error, since different variables have - % differnet dynamics and quantities - infcof = max(abs(Xh(:,end)./(x_new+eta))); - trc = infcof^(1/K); - LTE = infcof/(1-trc); -% LTE = infcof; -% LTE = max(abs(Xh(:,end)./x_new)); - LTE3 = max(abs(Xh(:,end-1)./(x_new+eta))); - - % relative error used to change order - LTE2 = infcof/min([norm(Xh(:,end-1)./Xh(:,end),Inf), sqrt(norm(Xh(:,end-2)./Xh(:,end),Inf)), sqrt(norm(Xh(:,end-3)./Xh(:,end-1),Inf))]); -% LTE2 = min(abs(Xh(:,end)./abs(Xh(:,end-1)))) -end \ No newline at end of file diff --git a/EMT/multiwindow.m b/EMT/multiwindow.m deleted file mode 100644 index a530355..0000000 --- a/EMT/multiwindow.m +++ /dev/null @@ -1,375 +0,0 @@ -function [t,x,y,hh,KK,bx,by,bz] = multiwindow(daeModelEquations,simu_T,x0,GLO,method,solverOptions,bx0) -format long -% multiwindow performs general multi-stage numerical integration. It -% applies to RK4, DT or other use-specified numerical solvers -% -% Inputs -% daeModelEquations the struct for equations related to a DAE model -% simu_T the simulation window: start, step length, end -% x0 initial values of the DAE -% GLO parameters in the DAE model -% method solution method -% solverOptions settings and parameters of the solution method -% Outputs -% t vector of time steps in the whole simulation -% x trajectory of state variables -% y trajectory of algebraic variables -% hh the time step used in each simulation step -% KK the order used in each simulation step -% -% See also Parareal, rk4, DT, getDAEModelEquations, getDAESolver - -f = daeModelEquations.fun_DAE; -fun_idx = daeModelEquations.fun_idx; - -%% -K = solverOptions.K; -% FS1_VS0 = solverOptions.FS1_VS0; -% FO1_VO0 = solverOptions.FO1_VO0; -hmax = solverOptions.hmax; -hmin = solverOptions.hmin; -Kmax = solverOptions.Kmax; -Kmin = solverOptions.Kmin; -% eps1 = solverOptions.eps1; -% eps2 = solverOptions.eps2; -% eps3 = solverOptions.eps3; -% eps4 = solverOptions.eps4; -tol = solverOptions.RBtol; -% q1 = solverOptions.q1; -% q2 = solverOptions.q2; -% dK = solverOptions.dK; -%% -a = simu_T(1); % start time of simulation -b = simu_T(2); % end time of simulation -h = simu_T(3); % time step length - -%% Judge if x0 is column vector or not. If rox vector, convert to column vector. -m = size(x0,1); % number of rows -if m == 1 - x0 = x0'; - m =size(x0,1); -end -bm = size(bx0,1); % number of rows -if bm == 1 - bx0 = bx0'; -end -%% -fun_y = str2func(strcat(func2str(f),'_AE')); -idx = fun_idx(GLO); -y0 = fun_y(x0,idx,GLO); -by0 = fun_y(bx0,idx,GLO); -%% Initial values of t and x -% mite = ceil((b-a)./hmin); -mite = 20000; -bt = zeros(1,mite); -bx = zeros(m,mite); -by = zeros(size(y0,1),mite); -t = zeros(1,mite); -x = zeros(m,mite); -y = zeros(size(y0,1),mite); -hh = zeros(1,mite); -KK = zeros(1,mite); -t(1) = a; -x(:,1) = x0; %initial conditions -bt(:,1) = a; -bx(:,1) = bx0; - -y(:,1) = y0; -by(:,1) = by0; -hh(:,1)= h; -KK(:,1)= K; - -%% Implement the integration method -ratio_b4= nan; -theta_max=2; - - - -ite = 1 ; - -% K = -(1/2)*log(tol); -M = 1; -fac1= 0.98; -fac2 = 1; - -% fac1 = 0; -% fac2 = 0; -% - -p= 1 ; -pf1 = 2*GLO.machine; -pf2 = 2*GLO.machine+4*GLO.machine*GLO.machine+idx.Nx*(idx.Nx+idx.Ny)+(idx.Nx); -pf3 = (idx.Nx); -% pf1 = 17*GLO.machine/2; -% pf2 = 27*GLO.machine/2+4*GLO.machine*GLO.machine+idx.Nx*(idx.Nx+idx.Ny)+(idx.Nx); -% pf2 = 0; -% pf3 = 0; -while t(ite)+h <= b - -switch solverOptions.solver - case 'DT' - - [t_new,x_new,y_new, LTE,LTE2,LTE3] = method(t(ite),x(:,ite),f,GLO,h,K); - ratio = LTE/h; - K_I = 0.3/K; - K_p = 0.4/K; - if isnan(ratio_b4) - else - theta = ((tol/ratio)^(K_I))*((ratio_b4/ratio)^(K_p)); - - if mod(ite,M)==0 -% Lin = ((K+p)/(K)) -% Lin = ((K+p)^2/(K)^2) - Lin = (pf1*(K+p)^2+pf2*(K+p)+pf3)/(pf1*(K)^2+pf2*(K)+pf3); -% Lde = ((K-p)/(K)) -% Lde = ((K-p)^2/(K)^2) - Lde = (pf1*(K-p)^2+pf2*(K-p)+pf3)/(pf1*(K)^2+pf2*(K)+pf3); - -% trc_in =(LTE2); -% trc_de= (LTE3); - - trc_in =(LTE2)/(1-LTE2^((K+p))); - trc_de= (LTE3)/(1-LTE3^((K-p))); - - tol_in = tol; - tol_de = tol; - - theta_in =1*(h*tol_in/trc_in)^(1/(K+p)); - theta_de =1*(h*tol_de/trc_de)^(1/(K-p)); - -% -% theta_in =1*(tol_in/trc_in)^(1/(K+p)); -% theta_de =1*(tol_de/trc_de)^(1/(K-p)); - -% % -% % theta_in =0.9*(h*tol/trc_in)^(1/(K+p)); -% % theta_de =0.9*(h*tol/trc_de)^(1/(K-p)); - - if theta> theta_max - h_can = theta_max*h; - else - h_can = theta*h; - end - - - if h_can > hmax - h_can = hmax; - - elseif h < hmin - h_can = hmin; - - end - - if theta_in> theta_max - theta_in = theta_max; - end - if theta_de> theta_max - theta_de = theta_max; - end - - - - if h*theta_in > hmax - h_inp = hmax; - - elseif h*theta_in < hmin - h_inp = hmin; - else - - h_inp = h*theta_in; - end - - - if h*theta_de > hmax - h_dep = hmax; - - elseif h*theta_de < hmin - h_dep = hmin; - else - h_dep = h*theta_de; - end - - - - - - - - - -% if min(hh(end-1:-1:end-M))>= h*theta&&Lin = h_can||KK(end-1)K) -% K=K-p; -% theta =theta_de -% -% end - - - if (Lde=K)))&&((h>= h_can||(h==hmax&&KK(1,ite+1)<=K))&&Lin =K)) - K=K-p; - theta =theta_de; - tol = tol_de; - - elseif (h>= h_can||(h_can==hmax&&KK(1,ite+1)<=K))&&Lin theta_max - h = theta_max*h; - else - h = theta*h; - end - - - end - ratio_b4 = ratio; - - - - case 'RK4' - [t_new,x_new,y_new] = method(t(ite),x(:,ite),f,GLO,h,K); - -end - - - -% bhx = bx(:,ite); -% for bht =t(ite):0.00031: t_new(end) -% if t_new(end)-bht>=0.00031 -% [bht_new,bhx_new,bhy_new] = rk4(bht,bhx,f,GLO,0.00031,8); -% bhx = bhx_new(:,end); -% else -% [bht_new,bhx_new,bhy_new] = rk4(bht,bhx,f,GLO,t_new(end)-bht,8); -% break -% % bht = t_new(end); -% end -% end - - if ite+1 > mite - mite = mite + 20000; -% bt(:,mite) = 0; -% bx(:,mite) = 0; -% by(:,mite) = 0; - t(:,mite) = 0; - x(:,mite) = 0; - y(:,mite) = 0; - hh(:,mite) = 0; - KK(:,mite) = 0; - end - - - - if h > hmax - h = hmax; - - elseif h < hmin - h = hmin; - - end - - if K > Kmax; K = Kmax; end - if K < Kmin; K = Kmin; end -% bt(:,ite+1) = bht_new(end); -% bx(:,ite+1) = bhx_new(:,end); -% by(:,ite+1) = bhy_new(:,end); - t(:,ite+1) = t_new; - x(:,ite+1) = x_new; - y(:,ite+1) = y_new; - hh(:,ite+1) = h; - KK(:,ite+1) = K; - -% r = tol^(1/(1+K)); - -% theta =(tol*h/r_norm)^(1/K); - - -% K; - ite = ite+1; - -% if FS1_VS0 == 0 -% if LTE < eps1; h = h*q1; end -% if LTE > eps2; h = h*q2; end -% if h > hmax; h = hmax; end -% if h < hmin; h = hmin; end -% end -% -% if FO1_VO0 == 0 -% if LTE < eps3 -% K = K-dK; -% end -% if LTE > eps4 -% K = K+dK; -% end -% if K > Kmax; K = Kmax; end -% if K < Kmin; K = Kmin; end -% end - if b-t(ite)0 - h = b-t(ite); - end - -end - -%% transpose: use .' to avoid conjugate transpose in case some variables are complex -t = t(:,1:ite).'; -x = x(:,1:ite).'; -y = y(:,1:ite).'; -hh = hh(:,1:ite).'; -KK = KK(:,1:ite).'; -% bx = bx(:,1:ite).'; -% by = by(:,1:ite).'; -% bz =max(abs(x-bx),[],2); -bx = x; -by = y; -bz =max(abs(x-bx),[],2); -% bz =max((x-bx)./(max(bx+10^(-19),[],1)),[],2); -end \ No newline at end of file