Derivative after parameter - system of ode

Question

\begin{align} \frac{dx_{1}}{dt}& = x_{2} \\ \frac{dx_{2}}{dt} &= - A x_{1} - B x_{2} + a\sin(Ct) \end{align} $x_{1}(0)=1$ and $x_{2}(0)=0$

I want to draw 2 plots in Matlab: $\frac{dx_{1}}{dA}$ and $\frac{dx_{2}}{dA}$

I found: $\frac{dx_{1}}{dA}$ = t and $\frac{dx_{2}}{dA}$ = $\int_{0}^{t} x_{2} e^t ds$

How can i solve $\frac{dx_{2}}{dA}$ ? Can somebody tell me how i can do it in matlab?

Answer 1

Whenever you have a parameter depending ODE $$ \dot x=f(t,x,q),~~~x(0,q)=x_0(q) $$ then the derivation of the solution for the parameters takes the form $$ \frac∂{∂t}\frac∂{∂q_k}x(t,q)=\frac∂{∂x}f(t,x,q)\frac∂{∂q_k}x(t,q)+\frac∂{∂q_k}f(t,x,q), \\~\\ \frac∂{∂q_k}x(0,q)=\frac∂{∂q_k}x_0(q), $$ or perhaps more readable, set $u_k=\frac∂{∂q_k}x$ and use "operator" notation, then $$ ∂_tu_k(t,q)=∂_xf(t,x,q)u_k(t,q)+∂_{q_k}f(t,x,q), \\~\\ u_k(0,q)=∂_{q_k}x_0(q). $$ Note that in general this equation contains the original solution $x$, so it increases the system. In a numerical solver, the new equations can be solved along the original equations with a correspondingly increased state space vector.

---

In the current case, set $u=\frac{∂x}{∂A}$, then \begin{align} \dot u_1&=u_2,&~~u_1(0)=0,\\ \dot u_2&=-x_1-Au_1-Bu_2& u_2(0)&=0. \end{align} This in no way has the stated solutions. In matlab you would formulate the extended system like

function xdot=ext_sys(t,x,a,A,B,C)
  u = x(3:4);
  dx = [ x(2); -A*x(1)-B*x(2)+a*sin(C*t) ]
  du = [ u(2); -x(1)-A*u(1)-B*u(2) ]
  xdot = [ dx ; du ]

---

(Update 2/27/2021) Debugging Kutta's 3/8 method in the comments with better (slightly) formatting

function [t, x] = RK38th(f, x0, t0, T, h)
  t = t0:h:T;
  nt = numel(t);
  nx = numel(x0);
  x = nan(nx, nt);
  x(:,1) = x0;
  for k = 1:nt-1
    k1 = h*f(t(k), x(:,k));
    k2 = h*f(t(k) + h/3 ,  x(:,k) + k1/3 );
    k3 = h*f(t(k) + (2*h)/3, x(:,k) - k1/3 + k2);
    k4 = h*f(t(k) + h, x(:,k) + k1 - k2 + k3);
    dx = (k1 + 3*k2 + 3*k3 + k4)/8;
    x(:,k+1) = x(:,k) + dx;
  end
end

Note that in k2, the x update is now like in the other stages composed of the k vectors.

This then results in the following plot, where the left part shows the solution, and the right part compares the computed derivative with a difference quotient, where the difference in $A$ was chosen large as $0.1$ to get visible deviations. For smaller values of eps, the curves are visually the same.

The parameters and calling code are

  A = 2; B = 0.02; C = 1; a = 0.8;
  t0 = 0; T = 10; h = 0.5;

  x0 = [2; 0];
  x0_ext = [ x0(1); x0(2); 0; 0 ];
  
  eps = 1e-1;
  [T0,X0] = RK38th(@(t,x)ext_sys(t,x,a,A,B,C), x0_ext, t0, T, h);
  [T1,X1] = RK38th(@(t,x)ext_sys(t,x,a,A+eps,B,C), x0_ext, t0, T, h);

  subplot(121);
  plot(T0,X0);
  legend(['x1'; 'x2'; 'u1'; 'u2']);

  subplot(122);
  plot(T0,(X1(1:2,:)-X0(1:2,:))/eps, '-+', T0, X0(3:4,:));
  legend(['dx1'; 'dx2'; 'u1'; 'u2']);