how to calculate the drivative of discretized ODE
1 visualización (últimos 30 días)
Mostrar comentarios más antiguos
Muhammad
el 6 de Dic. de 2023
Comentada: Muhammad
el 6 de Dic. de 2023
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)
0 comentarios
Respuesta aceptada
Torsten
el 6 de Dic. de 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 comentarios
Torsten
el 6 de Dic. de 2023
Editada: Torsten
el 6 de Dic. de 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
Más respuestas (0)
Ver también
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!