Dormand/Prince Runge-Kutta de orden 4 y 5

El Dopri54 es un método de Dormand/Prince Runge-Kutta cuyos pares están formados por Runge-Kutta de orden 4 y 5. Esto significa que se usa el método de orden 5 para avanzar cada paso y el método de orden 4 para estimar el error.

Elegimos el problema de valores iniciales

y'=y
y(0)=1
queremos saber el valor de y en t = 1. La solución obvia es

y(t)=et, por tanto y(1) = e

Codigo en Matlab para el método Dormand/Prince Dopri54: File dopri54.m



function [y,stad]=dopri54c(funcion,t0,t1,y0,tol,rpar)
% DOPRI54 solves system of ODEs with the Dormand and Prince 54 code.
% INPUT
% funcion - String containing name of user-supplied
% problem description.
% Call: yprime = fun(t,y,rpar) where funcion = 'fun'.
% t - Time (scalar).
% y - Solution column-vector.
% rpar - parameters possibly needed in buinding yprime
% yprime - Returned derivative column-vector; yprime(i)
% = dy(i)/dt.
% t0 - Initial value of t.
% t1 - Final value of t (must be t1 > t0).
% y0 - Initial value column-vector.
% tol - tolerance for the local error
% OUTPUT
% y - computed solution at y=t1.
% stad - vector containing number of steps, rejections and
% function calls
t=t0;
y=y0;
h=tol^(1/5)/4;
step=0;
nrej=0;
fcall=1;
a4=[44/45 -56/15 32/9]';
a5=[19372/6561 -25360/2187 64448/6561 -212/729]';
a6=[9017/3168 -355/33 46732/5247 49/176 -5103/18656]';
a7=[35/384 0 500/1113 125/192 -2187/6784 11/84]';
e=[71/57600 -1/40 -71/16695 71/1920 -17253/339200 22/525]';
k1=feval(funcion,t,y,rpar);

while t < t1
if t+h > t1; h=t1-t; end
k2=feval(funcion,t+h/5,y+h*k1/5,rpar);
k3=feval(funcion,t+3*h/10,y+h*(3*k1+9*k2)/40,rpar);
k4=feval(funcion,t+4*h/5,y+h*(a4(1)*k1+a4(2)*k2+a4(3)*k3),rpar);
k5=feval(funcion,t+8*h/9,y+h*(a5(1)*k1+a5(2)*k2+a5(3)*k3+...
a5(4)*k4),rpar);
k6=feval(funcion,t+h,y+h*(a6(1)*k1+a6(2)*k2+a6(3)*k3+a6(4)*k4+...
a6(5)*k5),rpar);
yt=y+h*(a7(1)*k1+a7(3)*k3+a7(4)*k4+a7(5)*k5+a7(6)*k6);
k2=feval(funcion,t+h,yt,rpar);
est=norm(h*(e(1)*k1+e(2)*k2+e(3)*k3+e(4)*k4+e(5)*k5+...
e(6)*k6),inf);
fcall=fcall+6;
% [t h est]

if est < tol
t=t+h;
k1=k2;
step=step+1;
y=yt;
else
nrej=nrej+1;
end
h=.9*min((tol/(est+eps))^(1/5),10)*h;

end
if nargout > 1
stad=[step nrej fcall];
end

Archivo funcion.m

Este es el código que contiene la función a evaluar, en este caso: y' = y

function[y] = funcion(x, y, rpar)
u=x;

Ejecución


Guardamos los archivos anteriores en el mismo directorio y desde la consola de Matlab nos movemos hasta él con el comando cd
Ahora ejecutamos desde la consola el comando
>>y = dopri54c('funcion', 0, 1, y0, 0.0001,0);

Ahora tenemos almacenado el resultado en la variable y, usamos
>> y
eso nos muestra el resultado obtenido

y = 2.718...

Ejectua aquí el método de Dormand/Prince 5_4 en The Runge kutta calculator





Ha sido util? Alguna idea para complementar el texto?



Deja tu post

Comentarios de otros usuarios





Deja tu post