how to modify code for system of delay differential equations with two delays
I have this code and i want to modify this for the system with multiple delays
%% ode simulation data, drop the transient part
tau = 1;
beta = 4;
n = 9.65;
gamma = 2;
trueModel = @(t,x,xdelay) beta*xdelay./(1+xdelay.^n)-gamma*x;
num_traj=100;
tr_c =linspace(0.5,1.5,num_traj); % 100 sets (ORIGINAL)
T = 20;
dt = 0.05;
tspan = [0 T];
N = round(T/dt)+1;
Tst = 10; % drop first 10 second
T_tr = 7;
N_tr = round(T_tr/dt)+1; % take fisrt part of remaining data for training
tau_max = 1.5; %%%FIXED DELAY FOR ODE DISCRETIZATION
M=30; %%%COLLOCATION DEGREE (ADDED), SHOULD CORRESPOND TO nx-1
[Dcheb,xcheb]=difmat(-tau_max,0,M);%%%CHEBYSHEV DIFFERENTIATION MATRIX AND NODES
wcheb=barywei(xcheb); %%%CHEBYSHEV BARYCENTIRC WEIGHTS
nX = round(tau_max/dt)+1; % number of states (OLD)
Xtrain_all = {};
BaseNet0 = struct;
BaseNet0.DiffMat = difmat(-tau_max,0,M);
BaseNet0.dt = dt;
% generate data + derivative
for k = 1:length(tr_c)
hist = @(t) tr_c(k);
sol = dde23(trueModel,tau,hist,tspan);
tint = linspace(0,T,N);
yint{k} = deval(sol,tint); %%%ALL SAMPLES IN [0,T]
tint_tr = linspace(Tst,Tst+T_tr,N_tr); %%%ALL TRAINING SAMPLES [10,17]
yint_tr{k} = yint{k}(round(Tst/dt):round(Tst/dt)+N_tr-1);
% obtain derivative via central difference
dy{k} = ctrDiff(yint{k},tint);
dy_tr{k} = ctrDiff(yint_tr{k},tint_tr);
m = length(tint_tr)-nX+1; % number of datapoints
xTrain = zeros(M+1,m);
dxTrain = zeros(M+1,m);
tchebTrain = zeros(M+1,m);
for i = 1:m
[~,tchebTrain(:,i)]=difmat(tint_tr(i),tint_tr(nX+i-1),M);%%%REINTERPOLATE AT CHEBYSHEV NODES
xTrain(:,i)=deval(sol,tchebTrain(:,i));
xTrainDelay(:,i)=deval(sol,tchebTrain(:,i)-tau);
dxTrain(:,i) = trueModel(1,xTrain(:,i),xTrainDelay(:,i));
end
Xtrain_all{k}=xTrain;
x0{k} = xTrain(:,1);
dXtrain_all{k}=dxTrain;
end
function v=barywei(x)
%BARYWEI barycentric weights.
% v=barywei(x) returns the weights v for barycentric interpolation at
% the nodes x according to [1].
%
% REFERENCES:
% [1] J.P. Berrut and L.N. Trefethen,"Barycentric Lagrange
% interpolation",SIAM Rev. 46(3):501-517,2004.
n=length(x); %number of nodes
dx=ones(n,1);
v=dx;
for m=2:n
for i=1:m-1
dx(i)=x(m)-x(i);
v(i)=-dx(i)*v(i);
end
v(m)=prod(dx(1:m));
end
v=1./v;
v=v/max(v); %normalization to 1
end
function [D,x]=difmat(a,b,N)
%DIFMAT Chebyshev nodes and pseudospectral differentiation matrix.
% [D,x]=DIFMAT(a,b,N) returns the pseudospectral differentiation
% matrix D of order N+1 on the N+1 Chebyshev nodes x in [a,b]
% according to [1].
% INPUT:
% a: left extremum (1×1)
% b: right extremum (1×1)
% N: degree of discretization (1×1)
% OUTPUT:
% D: differentiation matrix ((N+1)x(N+1))
% x: Chebyshev nodes ((N+1)x1)
%
% References:
% [1] L.N. Trefethen, "Spectral methods in Matlab", SIAM, 2000.
if N==0
x=1;
D=0;
return
end
x=(b+a+(b-a)*(cos(pi*(0:N)’/N)))/2;
c=[2;ones(N-1,1);2].*(-1).^(0:N)’;
X=repmat(x,1,N+1);
dX=X-X’;
D=(c*(1./c)’)./(dX+(eye(N+1)));
D=D-diag(sum(D’));
end
function dy = ctrDiff(yint,tint)
dt = tint(2)-tint(1);
dy = zeros(size(yint));
for k = 2:length(yint)-1
dy(:,k) = (yint(:,k+1)-yint(:,k-1))/2/dt;
end
dy(:,1)=(yint(:,2)-yint(:,1))/dt;
dy(:,end)=(yint(:,end)-yint(:,end-1))/dt;
end
my new model (system) is
tau = [1.5, 2];
par = [0.2, 0.5, 0.2,0.2,1];
trueModel =@(t,x,Z,par) [-x(3) – x(2)+ par(1) * Z(1,1)+ par(2) * Z(1,2);…
x(1) + par(3) * x(2);…
par(4) + x(3) * (x(1) – par(5))];
with history function
hist =@(t) kron(ones(length(t),1),[1.5; 0.4; 0.9]);
when i do modification for system i get an error in xTrain part
below is my error message
Unable to perform assignment because the size of the left side is 21-by-1 and the size of the right side is
3-by-21.
xTrain(:,i)=deval(sol,tchebTrain(:,i));I have this code and i want to modify this for the system with multiple delays
%% ode simulation data, drop the transient part
tau = 1;
beta = 4;
n = 9.65;
gamma = 2;
trueModel = @(t,x,xdelay) beta*xdelay./(1+xdelay.^n)-gamma*x;
num_traj=100;
tr_c =linspace(0.5,1.5,num_traj); % 100 sets (ORIGINAL)
T = 20;
dt = 0.05;
tspan = [0 T];
N = round(T/dt)+1;
Tst = 10; % drop first 10 second
T_tr = 7;
N_tr = round(T_tr/dt)+1; % take fisrt part of remaining data for training
tau_max = 1.5; %%%FIXED DELAY FOR ODE DISCRETIZATION
M=30; %%%COLLOCATION DEGREE (ADDED), SHOULD CORRESPOND TO nx-1
[Dcheb,xcheb]=difmat(-tau_max,0,M);%%%CHEBYSHEV DIFFERENTIATION MATRIX AND NODES
wcheb=barywei(xcheb); %%%CHEBYSHEV BARYCENTIRC WEIGHTS
nX = round(tau_max/dt)+1; % number of states (OLD)
Xtrain_all = {};
BaseNet0 = struct;
BaseNet0.DiffMat = difmat(-tau_max,0,M);
BaseNet0.dt = dt;
% generate data + derivative
for k = 1:length(tr_c)
hist = @(t) tr_c(k);
sol = dde23(trueModel,tau,hist,tspan);
tint = linspace(0,T,N);
yint{k} = deval(sol,tint); %%%ALL SAMPLES IN [0,T]
tint_tr = linspace(Tst,Tst+T_tr,N_tr); %%%ALL TRAINING SAMPLES [10,17]
yint_tr{k} = yint{k}(round(Tst/dt):round(Tst/dt)+N_tr-1);
% obtain derivative via central difference
dy{k} = ctrDiff(yint{k},tint);
dy_tr{k} = ctrDiff(yint_tr{k},tint_tr);
m = length(tint_tr)-nX+1; % number of datapoints
xTrain = zeros(M+1,m);
dxTrain = zeros(M+1,m);
tchebTrain = zeros(M+1,m);
for i = 1:m
[~,tchebTrain(:,i)]=difmat(tint_tr(i),tint_tr(nX+i-1),M);%%%REINTERPOLATE AT CHEBYSHEV NODES
xTrain(:,i)=deval(sol,tchebTrain(:,i));
xTrainDelay(:,i)=deval(sol,tchebTrain(:,i)-tau);
dxTrain(:,i) = trueModel(1,xTrain(:,i),xTrainDelay(:,i));
end
Xtrain_all{k}=xTrain;
x0{k} = xTrain(:,1);
dXtrain_all{k}=dxTrain;
end
function v=barywei(x)
%BARYWEI barycentric weights.
% v=barywei(x) returns the weights v for barycentric interpolation at
% the nodes x according to [1].
%
% REFERENCES:
% [1] J.P. Berrut and L.N. Trefethen,"Barycentric Lagrange
% interpolation",SIAM Rev. 46(3):501-517,2004.
n=length(x); %number of nodes
dx=ones(n,1);
v=dx;
for m=2:n
for i=1:m-1
dx(i)=x(m)-x(i);
v(i)=-dx(i)*v(i);
end
v(m)=prod(dx(1:m));
end
v=1./v;
v=v/max(v); %normalization to 1
end
function [D,x]=difmat(a,b,N)
%DIFMAT Chebyshev nodes and pseudospectral differentiation matrix.
% [D,x]=DIFMAT(a,b,N) returns the pseudospectral differentiation
% matrix D of order N+1 on the N+1 Chebyshev nodes x in [a,b]
% according to [1].
% INPUT:
% a: left extremum (1×1)
% b: right extremum (1×1)
% N: degree of discretization (1×1)
% OUTPUT:
% D: differentiation matrix ((N+1)x(N+1))
% x: Chebyshev nodes ((N+1)x1)
%
% References:
% [1] L.N. Trefethen, "Spectral methods in Matlab", SIAM, 2000.
if N==0
x=1;
D=0;
return
end
x=(b+a+(b-a)*(cos(pi*(0:N)’/N)))/2;
c=[2;ones(N-1,1);2].*(-1).^(0:N)’;
X=repmat(x,1,N+1);
dX=X-X’;
D=(c*(1./c)’)./(dX+(eye(N+1)));
D=D-diag(sum(D’));
end
function dy = ctrDiff(yint,tint)
dt = tint(2)-tint(1);
dy = zeros(size(yint));
for k = 2:length(yint)-1
dy(:,k) = (yint(:,k+1)-yint(:,k-1))/2/dt;
end
dy(:,1)=(yint(:,2)-yint(:,1))/dt;
dy(:,end)=(yint(:,end)-yint(:,end-1))/dt;
end
my new model (system) is
tau = [1.5, 2];
par = [0.2, 0.5, 0.2,0.2,1];
trueModel =@(t,x,Z,par) [-x(3) – x(2)+ par(1) * Z(1,1)+ par(2) * Z(1,2);…
x(1) + par(3) * x(2);…
par(4) + x(3) * (x(1) – par(5))];
with history function
hist =@(t) kron(ones(length(t),1),[1.5; 0.4; 0.9]);
when i do modification for system i get an error in xTrain part
below is my error message
Unable to perform assignment because the size of the left side is 21-by-1 and the size of the right side is
3-by-21.
xTrain(:,i)=deval(sol,tchebTrain(:,i)); I have this code and i want to modify this for the system with multiple delays
%% ode simulation data, drop the transient part
tau = 1;
beta = 4;
n = 9.65;
gamma = 2;
trueModel = @(t,x,xdelay) beta*xdelay./(1+xdelay.^n)-gamma*x;
num_traj=100;
tr_c =linspace(0.5,1.5,num_traj); % 100 sets (ORIGINAL)
T = 20;
dt = 0.05;
tspan = [0 T];
N = round(T/dt)+1;
Tst = 10; % drop first 10 second
T_tr = 7;
N_tr = round(T_tr/dt)+1; % take fisrt part of remaining data for training
tau_max = 1.5; %%%FIXED DELAY FOR ODE DISCRETIZATION
M=30; %%%COLLOCATION DEGREE (ADDED), SHOULD CORRESPOND TO nx-1
[Dcheb,xcheb]=difmat(-tau_max,0,M);%%%CHEBYSHEV DIFFERENTIATION MATRIX AND NODES
wcheb=barywei(xcheb); %%%CHEBYSHEV BARYCENTIRC WEIGHTS
nX = round(tau_max/dt)+1; % number of states (OLD)
Xtrain_all = {};
BaseNet0 = struct;
BaseNet0.DiffMat = difmat(-tau_max,0,M);
BaseNet0.dt = dt;
% generate data + derivative
for k = 1:length(tr_c)
hist = @(t) tr_c(k);
sol = dde23(trueModel,tau,hist,tspan);
tint = linspace(0,T,N);
yint{k} = deval(sol,tint); %%%ALL SAMPLES IN [0,T]
tint_tr = linspace(Tst,Tst+T_tr,N_tr); %%%ALL TRAINING SAMPLES [10,17]
yint_tr{k} = yint{k}(round(Tst/dt):round(Tst/dt)+N_tr-1);
% obtain derivative via central difference
dy{k} = ctrDiff(yint{k},tint);
dy_tr{k} = ctrDiff(yint_tr{k},tint_tr);
m = length(tint_tr)-nX+1; % number of datapoints
xTrain = zeros(M+1,m);
dxTrain = zeros(M+1,m);
tchebTrain = zeros(M+1,m);
for i = 1:m
[~,tchebTrain(:,i)]=difmat(tint_tr(i),tint_tr(nX+i-1),M);%%%REINTERPOLATE AT CHEBYSHEV NODES
xTrain(:,i)=deval(sol,tchebTrain(:,i));
xTrainDelay(:,i)=deval(sol,tchebTrain(:,i)-tau);
dxTrain(:,i) = trueModel(1,xTrain(:,i),xTrainDelay(:,i));
end
Xtrain_all{k}=xTrain;
x0{k} = xTrain(:,1);
dXtrain_all{k}=dxTrain;
end
function v=barywei(x)
%BARYWEI barycentric weights.
% v=barywei(x) returns the weights v for barycentric interpolation at
% the nodes x according to [1].
%
% REFERENCES:
% [1] J.P. Berrut and L.N. Trefethen,"Barycentric Lagrange
% interpolation",SIAM Rev. 46(3):501-517,2004.
n=length(x); %number of nodes
dx=ones(n,1);
v=dx;
for m=2:n
for i=1:m-1
dx(i)=x(m)-x(i);
v(i)=-dx(i)*v(i);
end
v(m)=prod(dx(1:m));
end
v=1./v;
v=v/max(v); %normalization to 1
end
function [D,x]=difmat(a,b,N)
%DIFMAT Chebyshev nodes and pseudospectral differentiation matrix.
% [D,x]=DIFMAT(a,b,N) returns the pseudospectral differentiation
% matrix D of order N+1 on the N+1 Chebyshev nodes x in [a,b]
% according to [1].
% INPUT:
% a: left extremum (1×1)
% b: right extremum (1×1)
% N: degree of discretization (1×1)
% OUTPUT:
% D: differentiation matrix ((N+1)x(N+1))
% x: Chebyshev nodes ((N+1)x1)
%
% References:
% [1] L.N. Trefethen, "Spectral methods in Matlab", SIAM, 2000.
if N==0
x=1;
D=0;
return
end
x=(b+a+(b-a)*(cos(pi*(0:N)’/N)))/2;
c=[2;ones(N-1,1);2].*(-1).^(0:N)’;
X=repmat(x,1,N+1);
dX=X-X’;
D=(c*(1./c)’)./(dX+(eye(N+1)));
D=D-diag(sum(D’));
end
function dy = ctrDiff(yint,tint)
dt = tint(2)-tint(1);
dy = zeros(size(yint));
for k = 2:length(yint)-1
dy(:,k) = (yint(:,k+1)-yint(:,k-1))/2/dt;
end
dy(:,1)=(yint(:,2)-yint(:,1))/dt;
dy(:,end)=(yint(:,end)-yint(:,end-1))/dt;
end
my new model (system) is
tau = [1.5, 2];
par = [0.2, 0.5, 0.2,0.2,1];
trueModel =@(t,x,Z,par) [-x(3) – x(2)+ par(1) * Z(1,1)+ par(2) * Z(1,2);…
x(1) + par(3) * x(2);…
par(4) + x(3) * (x(1) – par(5))];
with history function
hist =@(t) kron(ones(length(t),1),[1.5; 0.4; 0.9]);
when i do modification for system i get an error in xTrain part
below is my error message
Unable to perform assignment because the size of the left side is 21-by-1 and the size of the right side is
3-by-21.
xTrain(:,i)=deval(sol,tchebTrain(:,i)); differential equations, delay differential equations MATLAB Answers — New Questions