Prof. Vladimir A. Dobrushkin
This tutorial contains many matlab scripts.
You, as the user, are free to use all codes for your needs, and have the right to distribute this tutorial and refer to this tutorial as long as this tutorial is accredited appropriately. Any comments and/or contributions for this tutorial are welcome; you can send your remarks to
The derivatives y' and y'' must be extressed in terms of f(x,y) and its partial derivatives according to the given differential equation \( y' = f(x,y) . \) So we have
Since there are only three equations in four unknowns, the system is underdetermined and we are permitted to choose one of the coefficients as we wish. There are several special choices that have been studied in the literature, and we mention some of them.
If we choose b1 = 1/2, this yields b2 = 1/2, c2 = 1, and a2,1 = 1. So we get Heun's method.
If we choose b1 = 0, this yields b2 = 1, c2 = 1/2, and a2,1 = 1/2. So we get the modified Euler method.
Example. The Butcher table for the explicit midpoint method is given by:
function Complete
Complete
I. Second Order Optimal Runge-Kutta Approximation
Ralston's method is a second-order method with two stages and a minimum local error bound. Its Mathematica realization is presented below when the step size is denoted by h: \begin{equation} \label{E35.opt} y_{n+1} =y_n + \frac{h}{4}\,f(x_n , y_n ) + \frac{3h}{4} \,f\left( x_n + \frac{2h}{3} , y_n + \frac{2h}{3} \,f(x_n , y_n ) \right) , \quad n=0,1,2,\ldots . \end{equation} Its Butcher tableau is
II. Implicit Runge-Kutta of order 2
\begin{equation} \label{E36.implicit} y_{n+1} = y_n + \frac{h}{2}\left( k_1 + k_2 \right) , \end{equation} where \begin{align*} k_1 &= f\left( x_n + \frac{\sqrt{3} -1}{2\sqrt{3}}\,h , \ y_n + \frac{h}{4}\,k_1 + \frac{\sqrt{3} -2}{4\sqrt{3}}\,hk_2\right) , \\ k_2 &= f\left( x_n + \frac{\sqrt{3} +1}{2\sqrt{3}}\,h , \ y_n + \frac{\sqrt{3} +2}{4\sqrt{3}}\,hk_1 + \frac{h}{4}\,k_2\right) , \end{align*}
III. Ralston's method
The method was proposed by Anthony Ralston (born in 1930, in New York City). Anthony Ralston received his Ph.D. from the Massachusetts Institute of Technology in 1956 under the direction of M. Eric Reissner. For approximating the solution of the initial value problem
The second order Ralson approximation can written as follows:
\begin{equation*} %\label{E35.ralson} y_{n+1} =y_n + \frac{h}{3}\,f(x_n , y_n ) + \frac{2h}{3} \,f\left( x_n + \frac{3h}{4} , y_n + \frac{3h}{4} \,f(x_n , y_n ) \right) , \quad n=0,1,2,\ldots . \end{equation*}
IV. Generic second-order method
There is a family of Runge--Kutta methods of order 2, parameterized by an arbitrary positive number α from the interval (0,1] and given by the formula
function [y, t] = implicit_rk2_ode(f, t0, y0, h, N)
% RK4_ODE Solve a first order diff eq with implicit second order Runge-Kutta method
% [y, t] = IMPLICIT_RK4_ODE(f, t0, y0, h, N) solves the IVP
% {y' = f(t,y), y(t0) = y0} on the interval [t, t + h*N]
% Example usage:
% We want to solve the IVP: y' = 1 - t + 4y, y(0) = 1
% f = @(t,y) 1 - t + 4*y;
% t0 = 0; y0 = 1; h = 0.001; N = 1000;
% [y, t] = implicit_rk2_ode(f, t0, y0, h, N);
% plot(t,y)
% format long g
% disp [t' y']
t = t0:h:(t0+N*h); % set up t-points
y(N+1)=0; % memory allocation for y-points
% Important note: because Matlab starts with index 1, what we refer to as
% t_j is called t(j +1) in Matlab; similarly, what we refer to as y_j
% is called y(j +1) in Matlab
% First y-point is from our initial condition:
y(1) = y0;
% Now, the main loop of the algorithm:
for j = 0:N-1
t_j = t(j +1);
y_j = y(j +1);
% First, find k1 and k2:
[k1, k2] = solve_implicit_system(f, h, t_j, y_j);
% Now use the Runge-Kutta equation to find y_j+1:
y_jplus1 = y_j + (h/2)*(k1+k2);
y(j+2) = y_jplus1;
end
end
f = @(t,y) 1 - t + 4*y;
t0 = 0; y0 = 1; h = 0.001; N = 1000;
[y, t] = implicit_rk2_ode(f, t0, y0, h, N);
plot(t,y)
format long g
disp([t' y'])
y(t==1.0)
function rukun
f=inline('-x^2+y') % the right side of ODE
h=0.1;
x0=0; % intial argument
y0=-1; %initial value y(x0)
xn=1; % terminal argument
x=x0:h:xn;
y(numel(x))=0; %memory allocation
y(1)=y0;
for i=1:numel(x)-1
k=ss(f, h, x(i), y(i));
y(i+1)=y(i)+1/2*h*(k(1)+k(2));
end
disp([x',y'])
plot(x,y)
function solution=ss(f,h,x,y)
s=sqrt(3);
t1=(s-1)/(2*s);
t2=(s-2)/(4*s);
t3=(s+1)/(2*s);
t4=(s+2)/(4*s);
F=@(k) [f(x+t1*h, y+(h/4)*k(1)+t2*h*k(2))-k(1);
f(x+t3*h, y+t4*h*k(1)+h/4*k(2))-k(2)];
initial_guess=[1;1];
[solution]=fsolve(F, initial_guess);