Clear Filters
Clear Filters

how to calculate the drivative of discretized ODE

4 views (last 30 days)
I am discretizing DDE to ODE using pseudospectral method. I want to compute derivative of its solution for training state and want to use the right hand equation of the discretized ODE here dODE represents discretized ODE
tspan=[0 20]; M=5; t_r=0.8;
phi = @(x) cos(x);
g = @(t,y,Z,par) par * y * (1 - Z);
tau = 1;
par = 1.5;
[D, theta] = difmat(-tau, 0, M);
X = phi(theta);
u0=X;
options = odeset('RelTol', 1e-10, 'AbsTol', 1e-10);
sol = ode45(@(t,u) dODE(t, u, par, tau, M,D), tspan, u0,options);
x_n = sol.y ;
t = linspace(tspan(1), tspan(2), 100);
x_an = interp1(sol.x, x_n(1,:), t, 'linear');
n_s_r = size(x_an, 2);
x_a = floor(n_s_r * t_r);
x_tn = x_an(:, 1:x_a);
n_all = length(t);
n_train = round(t_r * n_all);
t_t = t(1:n_train);
function dydt = dODE(t,u, par, tau, M, D)
yM = u(1);
VM = u(2:end);
dMDM_DDE = kron(D(2:end,:), eye(1));
dydt = [par*yM*(1-VM(end)); (dMDM_DDE)*[yM;VM]];
end
function [D, x] = difmat(a, b, M)
% CHEB pseudospectral differentiation matrix on Chebyshev nodes.
% [D,x]=CHEB(a,b,M) returns the pseudospectral differentiation matrix D
if M == 0
x = 1;
D = 0;
return
end
x = ((b - a) * cos(pi * (0:M)' / M) + b + a) / 2;
c = [2; ones(M-1, 1); 2].*(-1).^(0:M)';
X = repmat(x, 1, M+1);
dX = X - X';
D = (c * (1./c)')./(dX + (eye(M+1)));
D = D - diag(sum(D'));
end
(like this way
for j = 1:length(t_t)
DX(:,j) = g(t_t(j), x_tn(j), x_dn(:,j), par);
end this is derivative of original DDE)

Accepted Answer

Torsten
Torsten on 6 Dec 2023
After the line
sol = ode45(@(t,u) dODE(t, u, par, tau, M,D), tspan, u0,options);
you can compute the derivatives of the solution as
[~,yp] = deval(sol,sol.x)
  3 Comments
Torsten
Torsten on 6 Dec 2023
Edited: Torsten on 6 Dec 2023
Here is an example:
fun = @(t,y) [y(1);-y(2)];
tspan = 0:0.1:1;
y0 = [1;1];
sol = ode45(fun,tspan,y0);
% Compute 1st derivative of the solution
[~,yp] = deval(sol,sol.x);
% Compute 2nd derivatve of the solution
n = size(yp,2);
ypp(:,1)=(yp(:,2) - yp(:,1))./(sol.x(2)-sol.x(1));
ypp(:,2:n-1) = (yp(:,3:n) - yp(:,1:n-2))./(sol.x(3:n)-sol.x(1:n-2));
ypp(:,n)=(yp(:,n) - yp(:,n-1))/(sol.x(n) - sol.x(n-1));
% Plot the 2nd derivative of the solution
plot(sol.x,ypp)
grid on

Sign in to comment.

More Answers (0)

Categories

Find more on Programming in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!