Implementing implicit method for system of ODEs

1 view (last 30 days)
Romio
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)')

Answers (1)

Torsten
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
Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient.
plot(t,y)
xlabel('t')
ylabel('y(t)')

Products


Release

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by