how to calculate the drivative of discretized ODE

1 次查看(过去 30 天)
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)

采纳的回答

Torsten
Torsten 2023-12-6
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 个评论
Torsten
Torsten 2023-12-6
编辑:Torsten 2023-12-6
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

请先登录,再进行评论。

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Programming 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by