# Implementing implicit method for system of ODEs

1 view (last 30 days)
Romio on 25 Jun 2022
Edited: Torsten on 26 Jun 2022
I am trying to implement implicit Runge Kutta for solving a system of ODEs. Currently I was able to do it for scalar ODE (y'=-y), but I got stuck on how to modify it so that it work for a system of ODEs. Something like y1' = y2, y2' = -y1.
Everythin before the loop is fine if we replace the function handle f by the following function
function dydt = f(t,y)
dydt(1) = y(2);
dydt(2) = -y(1);
end
and modify the initial value to be something like
y0 = [1;0]
But How do we modify the process of solving the nonlinear system inside the loop? Also I feel I have hard-coded the method for this specific 2 stage method. Is there any way to generalize the loop for any implicit runge kutta method with size s by s?
Here is the code
n = 50; % number of time steps
t0 = 0; % initial time
T = 10; % final time
dt = (T-t0)/n; % size of timestep
y0 = 1; % initial data
m = length(y0); % dim. of problem
% intitalize time and solution vector
t = zeros(n+1,1);
y = zeros(n+1,m);
% initial conditions
t(1,1) = t0;
y(1,:) = y0(:);
% RHS of ODE
f = @(t,y) -1 * y;
% Runge Kutta coefficients
A = [1/2 + 1/(2*sqrt(3)) , 0;
-1/sqrt(3), 1/2 + 1/(2*sqrt(3))];
b = [1/2;1/2];
c = [1/2 + 1/(2*sqrt(3));1/2 - 1/(2*sqrt(3))];
s = length(c);
for i = 1:n
XI = @(xi) [xi(1) - (y(i,:) + dt*A(1,1)*f(t(i)+c(1)*dt,xi(1))...
+ dt*A(1,2)*f(t(i)+c(2)*dt,xi(2)));
xi(2) - (y(i,:) + dt*A(2,1)*f(t(i)+c(1)*dt,xi(1))...
+ dt*A(2,2)*f(t(i)+c(2)*dt,xi(2)))];
xi = fsolve(XI,zeros(s,1));
t(i+1,1) = t(i,1) + dt;
y(i+1,:) = y(i,:) + dt * (b(1)*f(t(i)+c(1)*dt,xi(1))+...
b(2)*f(t(i)+c(2)*dt,xi(2)));
end
plot(t,y)
xlabel('t')
ylabel('y(t)')

Torsten on 26 Jun 2022
Edited: Torsten on 26 Jun 2022
n = 50; % number of time steps
t0 = 0; % initial time
T = 2*pi; % final time
dt = (T-t0)/n; % size of timestep
y0 = [1;0]; % initial data
m = length(y0); % dim. of problem
% intitalize time and solution vector
t = zeros(1,n+1);
y = zeros(m,n+1);
% initial conditions
t(1) = t0;
y(:,1) = y0;
% RHS of ODE
f = @(t,y) [y(2);-y(1)];
% Runge Kutta coefficients
A = [1/2 + 1/(2*sqrt(3)) , 0;
-1/sqrt(3), 1/2 + 1/(2*sqrt(3))];
b = [1/2;1/2];
c = [1/2 + 1/(2*sqrt(3));1/2 - 1/(2*sqrt(3))];
s = length(c);
xiguess = zeros(m*s,1);
for i = 1:n
XI = @(xi) [xi(1:m) - f(t(i)+c(1)*dt,y(:,i)+dt*(A(1,1)*xi(1:m)+A(1,2)*xi(m+1:m*s)));
xi(m+1:m*s) - f(t(i)+c(2)*dt,y(:,i)+dt*(A(2,1)*xi(1:m)+A(2,2)*xi(m+1:m*s)))];
xi = fsolve(XI,xiguess);
t(i+1) = t(i) + dt;
y(:,i+1) = y(:,i) + dt * (b(1)*xi(1:m) + b(2)*xi(m+1:m*s));
xiguess = xi;
end
plot(t,y)
xlabel('t')
ylabel('y(t)') R2021b

### Community Treasure Hunt

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

Start Hunting!