Theta Method for system of ODEs

Question

So here is a question about a little Matlab code writing because I've never seen "theta method"...

Consider the $θ$-method $$y_{n+1} = y_n + θhf_n + (1 − θ)hf_{n+1}$$ where $0 ≤ θ ≤ 1$ and $f_n = f(t_n, y_n)$. Write a MATLAB code to implement the theta method for systems of ODEs. For $θ = 0$, $0.5$, $1$, use your code for solving \begin{aligned} y_1' &= −y_1 \\ y_2' &= −100 \left(y_2 − \sin(t)\right) + \cos(t) \end{aligned} for $0 ≤ t ≤ 1$, with initial value $y_1 = 1$, $y_2 = 2$. Try this for stepsizes $h = .01$ and $h = .05$.

Now I know how to do the Euler Code....

h = 0.05; %step size
y1 = 1; %IC Conditions Y1
x1 = 2; %IC Condition Y2
t0 = 0; %Initial time
t1 = 1; %Final time

nsteps = (t1-t0)/h; %number of steps

y(1) = y1; %Array for euler method Y1, -I.C Y1
x(1) = x1; %Array for euler method Y2  -I.C Y2
t(1) = t0 ; %array for euler method t

%% For loop and implementation of Euler Method

for n = 1:nsteps
   t(n+1) = t(n)+h; %time stepping euler method, increment time
   yprime(n+1) = -y(n); %Right hand side of diff eqy y1' = -y1
   xprime(n+1) = -100*(x(n)-sin(t(n)))+cos(t(n)); %RHS of y2'= -100*(y2-sin(t))+cos(t)
   y(n+1) = y(n)+yprime(n+1)*h;   
   x(n+1) = x(n)+xprime(n+1)*h;
end

Would I just be changing the last part of the code where I have this?

   y(n+1) = y(n)+yprime(n+1)*h;   
   x(n+1) = x(n)+xprime(n+1)*h;

Essentially I would think I change that part above into the theta method of,

$y_{n+1} = y_n + θhf_n + (1 − θ)hf_{n+1}$

If anyone is familiar with this method, please let me know how you would write this code, any feedback is helpful.

Thank you

Answer 1

Let us introduce the vector of unknowns $Y = (y_1,y_2)^\top$ which satisfies the ODE system $Y'(t) = F(t,Y)$. At each step, one needs to solve $$ Y_{n+1} = Y_{n} + h\left(\theta\, F(t_n,Y_n) + \left(1-\theta\right) F(t_{n+1},Y_{n+1}) \right) , $$ where $t_n = nh$ and $Y_n \simeq Y(t_n)$. In general, this is a nonlinear equation for $Y_{n+1}$, which is rewritten $G(Y_{n+1}) = 0$, where $$ G : Y \mapsto Y - h\left( \left(1-\theta\right) F(nh+h,Y) + \theta\, F(nh,Y_n) \right) - Y_{n} \, . $$ It is solved at each step using Newton's method $$ Y_{n+1}^{k+1} = Y_{n+1}^{k} - \left( G'(Y_{n+1}^{k})\right)^{-1}\, G(Y_{n+1}^{k}) \, , $$ where $$ G' : Y \mapsto I - \left(1-\theta\right) h\, \partial_Y F(nh+h,Y) $$ is the Jacobian matrix of $G$, and $\partial_Y F(nh+h,\cdot)$ is the Jacobian matrix of $F(nh+h,\cdot)$. The initial guess is $Y_{n+1}^{0} = Y_n$, and the final value is $Y_{n+1}^\infty = Y_{n+1}$.

The code below performs the proposed integration:

% numerical parameters
theta = 0.5;
h = 0.01;
T = 1;

% functions describing ODE system
dim = 2;
F = @(t,Y) [-Y(1); -100*(Y(2)-sin(t))+cos(t)];
Fp = @(t,Y) [-1, 0; 0, -100];

% functions for Euler method
G = @(Y,Yn,n) Y - h*( (1-theta)*F(n*h+h,Y) + theta*F(n*h,Yn) ) - Yn;
Gp = @(Y,n) eye(dim) - h*(1-theta)*Fp(n*h+h,Y);

% initialization
N = floor(T/h);
Y = zeros(2,N);
Y(:,1) = [1; 2];
k = zeros(1,N-1);

for n=1:N-1
    Yn = Y(:,n);
    Yk = Yn;
    err = 1;
    while (err>1e-14)&&(k(n)<20)
        Ykp1 = Yk - Gp(Yk,n)\G(Yk,Yn,n);
        err = norm(Ykp1-Yk);
        Yk = Ykp1;
        k(n) = k(n)+1;
    end
    Y(:,n+1) = Ykp1;
end

plot(Y(1,:),Y(2,:),'bo');

% analytical solution
hold on
t = (0:N)*h;
plot(exp(-t),2*exp(-100*t)+sin(t),'k');

A comparison of the numerical solution (blue) with the analytical solution (red) in the phase space $y_1$-$y_2$ is given in the graphics output