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)
댓글 수: 0
채택된 답변
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
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개)
참고 항목
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!