| |
| Autor |
 2(3)-Runge-Kutta-Verfahren |
|
Specialagent
Aktiv 
Dabei seit: 17.11.2012
Mitteilungen: 116
 |  Themenstart: 2020-03-02
|
Hallo,
ich habe versucht das 2(3)-RK Verfahren zu implementieren. Gelöst werden soll die Prothero-Robinson DGL: y' = \lambda(y-g(t)) + g'(t). Dabei sei g(t) = exp(-t) und y_0 = g(0).
Ich soll das ganze für \lambda = -5, \lambda = -500 sowie für verschiedene Schrittweiten h = 2^(-k), k = 4,5,...,10 durchspielen.
Für \lambda = -5 läuft alles wie es zu erwarten war, aber leider nicht für \lambda = -500. Dort läuft das Programm ewig und ich muss es abbrechen.
Sieht jemand meinen Fehler? Vielen Dank für jede Hilfe.
P.S. Ich habe es mal ins Matlab Forum gepostet. Evtl. passt das aber auch ins Numerik Forum?
\sourceon matlab
function [x, y] = rk23(lambda,k,a,b,y0,rtol)
%Ausgabe:
%1) Diskretisierungspunkte x und numerisch bestimmten Werte y
%2) Grafische Ausgabe der exakten und numerischen Lösung
%Übergabeparameter:
%1) lambda aus der Prothero-Robinson DGL
%2) k für die Wahl der Schrittweite h = 2^-k, k = 4,5,...,10.
%3) Intervall [a,b]
%4) Anfangswert y0 der DGL
%5) relative Toleranz rtol jedes Schrittes für die Schrittweitensteuerung
f = @(t,y) lambda*(y - exp(-t)) - exp(-t);
h0 = 2^-k;
atol = 1e-13; %absolute Toleranz
alpha = 0.8; % 0 < alpha < 1 für Schrittweitensteuerung
% Anfangszeit
i = 1;
x(1) = a;
t = a;
% Anfangsbedingung
y(1,:) = y0;
y_n = y0;
% Wenn letzte Interation erreicht wird, setze iteration = 1, sonst = 0
iteration = 0;
while iteration == 0 %solange ende des Intervalls nicht erreicht ist,...
%Überprüfe ob Ende des Intervalls erreicht
if t + h0 > b
h0 = b - t;
iteration = 1;
end
% Berechnung der benötigten k-Werte
k1 = f(t,y_n);
k2 = f(t + 0.5*h0, y_n + 0.5*h0*k1);
k3 = f(t + h0, y_n + - h0*k1 + 2*h0*k2);
w = y_n + h0 * k2; %Verfahren von Collatz (Ordnung 2)
z = y_n + 1/6 * h0 * (k1 + 4*k2 + k3); %Verfahren von Kutta (Ordnung 3)
lf = w - z; %Fehler
% Prüfung der Toleranz
Tol = rtol * norm(y_n) + atol;
if abs(lf) <= Tol %abs(lf) <= rtol ?
t = t + h0;
h0 = nthroot(rtol/lf, 3)*h0; %Neue Schrittweite
i = i + 1;
x(i) = t;
y_n = z; %Speichere Werte von Ordnung 3
y(i,:) = z;
else
h0 = alpha * h0;
iteration = 0;
end
end
%Exakte Lösung
y_e = a:0.01:b;
y_exakt = exp(-y_e);
figure(3);
plot(y_e,y_exakt,'b','LineWidth',2);
hold on
plot(x,y,'go','MarkerSize',3,'LineWidth',2);
title('2(3) RK-Verfahren');
legend('Exakte Lsg', '2(3) RK');
\sourceoff
|
 Profil
|
Delastelle
Senior 
Dabei seit: 17.11.2006
Mitteilungen: 2320
 | Beitrag No.1, eingetragen 2020-03-03
|
Hallo Specialagent!
Eventuell kannst Du Dir mal die gemachten Schritte anschauen
(Werte bei einzelnen Schritten ausgeben).
Irgendwo hatte ich auch mal den Code von einem älteren Matlab ode23.m.
Ich muss mal suchen...
Teste eventuell auch mal Dein Programm für eine größere Toleranz (nicht 1e-13 sondern etwa 1e-6).
Viele Grüße
Ronald
Edit:
das Original-Ode23.m von Octave
\showon ode23.m Octave
%# Copyright (C) 2006-2012, Thomas Treichl
%# OdePkg - A package for solving ordinary differential equations and more
%#
%# This program is free software; you can redistribute it and/or modify
%# it under the terms of the GNU General Public License as published by
%# the Free Software Foundation; either version 2 of the License, or
%# (at your option) any later version.
%#
%# This program is distributed in the hope that it will be useful,
%# but WITHOUT ANY WARRANTY; without even the implied warranty of
%# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
%# GNU General Public License for more details.
%#
%# You should have received a copy of the GNU General Public License
%# along with this program; If not, see .
%# -*- texinfo -*-
%# @deftypefn {Function File} {[@var{}] =} ode23 (@var{@@fun}, @var{slot}, @var{init}, [@var{opt}], [@var{par1}, @var{par2}, @dots{}])
%# @deftypefnx {Command} {[@var{sol}] =} ode23 (@var{@@fun}, @var{slot}, @var{init}, [@var{opt}], [@var{par1}, @var{par2}, @dots{}])
%# @deftypefnx {Command} {[@var{t}, @var{y}, [@var{xe}, @var{ye}, @var{ie}]] =} ode23 (@var{@@fun}, @var{slot}, @var{init}, [@var{opt}], [@var{par1}, @var{par2}, @dots{}])
%#
%# This function file can be used to solve a set of non--stiff ordinary differential equations (non--stiff ODEs) or non--stiff differential algebraic equations (non--stiff DAEs) with the well known explicit Runge--Kutta method of order (2,3).
%#
%# If this function is called with no return argument then plot the solution over time in a figure window while solving the set of ODEs that are defined in a function and specified by the function handle @var{@@fun}. The second input argument @var{slot} is a double vector that defines the time slot, @var{init} is a double vector that defines the initial values of the states, @var{opt} can optionally be a structure array that keeps the options created with the command @command{odeset} and @var{par1}, @var{par2}, @dots{} can optionally be other input arguments of any type that have to be passed to the function defined by @var{@@fun}.
%#
%# If this function is called with one return argument then return the solution @var{sol} of type structure array after solving the set of ODEs. The solution @var{sol} has the fields @var{x} of type double column vector for the steps chosen by the solver, @var{y} of type double column vector for the solutions at each time step of @var{x}, @var{solver} of type string for the solver name and optionally the extended time stamp information @var{xe}, the extended solution information @var{ye} and the extended index information @var{ie} all of type double column vector that keep the informations of the event function if an event function handle is set in the option argument @var{opt}.
%#
%# If this function is called with more than one return argument then return the time stamps @var{t}, the solution values @var{y} and optionally the extended time stamp information @var{xe}, the extended solution information @var{ye} and the extended index information @var{ie} all of type double column vector.
%#
%# For example, solve an anonymous implementation of the Van der Pol equation
%#
%# @example
%# fvdb = @@(vt,vy) [vy(2); (1 - vy(1)^2) * vy(2) - vy(1)];
%#
%# vopt = odeset ("RelTol", 1e-3, "AbsTol", 1e-3, \
%# "NormControl", "on", "OutputFcn", @@odeplot);
%# ode23 (fvdb, [0 20], [2 0], vopt);
%# @end example
%# @end deftypefn
%#
%# @seealso{odepkg}
%# ChangeLog:
%# 20010703 the function file "ode23.m" was written by Marc Compere
%# under the GPL for the use with this software. This function has been
%# taken as a base for the following implementation.
%# 20060810, Thomas Treichl
%# This function was adapted to the new syntax that is used by the
%# new OdePkg for Octave and is compatible to Matlab's ode23.
function [varargout] = ode23 (vfun, vslot, vinit, varargin)
if (nargin == 0) %# Check number and types of all input arguments
help ('ode23');
error ('OdePkg:InvalidArgument', ...
'Number of input arguments must be greater than zero');
elseif (nargin < 3)
print_usage;
elseif ~(isa (vfun, 'function_handle') || isa (vfun, 'inline'))
error ('OdePkg:InvalidArgument', ...
'First input argument must be a valid function handle');
elseif (~isvector (vslot) || length (vslot) < 2)
error ('OdePkg:InvalidArgument', ...
'Second input argument must be a valid vector');
elseif (~isvector (vinit) || ~isnumeric (vinit))
error ('OdePkg:InvalidArgument', ...
'Third input argument must be a valid numerical value');
elseif (nargin >= 4)
if (~isstruct (varargin{1}))
%# varargin{1:len} are parameters for vfun
vodeoptions = odeset;
vfunarguments = varargin;
elseif (length (varargin) > 1)
%# varargin{1} is an OdePkg options structure vopt
vodeoptions = odepkg_structure_check (varargin{1}, 'ode23');
vfunarguments = {varargin{2:length(varargin)}};
else %# if (isstruct (varargin{1}))
vodeoptions = odepkg_structure_check (varargin{1}, 'ode23');
vfunarguments = {};
end
else %# if (nargin == 3)
vodeoptions = odeset;
vfunarguments = {};
end
%# Start preprocessing, have a look which options are set in
%# vodeoptions, check if an invalid or unused option is set
vslot = vslot(:).'; %# Create a row vector
vinit = vinit(:).'; %# Create a row vector
if (length (vslot) > 2) %# Step size checking
vstepsizefixed = true;
else
vstepsizefixed = false;
end
%# Get the default options that can be set with 'odeset' temporarily
vodetemp = odeset;
%# Implementation of the option RelTol has been finished. This option
%# can be set by the user to another value than default value.
if (isempty (vodeoptions.RelTol) && ~vstepsizefixed)
vodeoptions.RelTol = 1e-6;
warning ('OdePkg:InvalidArgument', ...
'Option "RelTol" not set, new value %f is used', vodeoptions.RelTol);
elseif (~isempty (vodeoptions.RelTol) && vstepsizefixed)
warning ('OdePkg:InvalidArgument', ...
'Option "RelTol" will be ignored if fixed time stamps are given');
end
%# Implementation of the option AbsTol has been finished. This option
%# can be set by the user to another value than default value.
if (isempty (vodeoptions.AbsTol) && ~vstepsizefixed)
vodeoptions.AbsTol = 1e-6;
warning ('OdePkg:InvalidArgument', ...
'Option "AbsTol" not set, new value %f is used', vodeoptions.AbsTol);
elseif (~isempty (vodeoptions.AbsTol) && vstepsizefixed)
warning ('OdePkg:InvalidArgument', ...
'Option "AbsTol" will be ignored if fixed time stamps are given');
else
vodeoptions.AbsTol = vodeoptions.AbsTol(:); %# Create column vector
end
%# Implementation of the option NormControl has been finished. This
%# option can be set by the user to another value than default value.
if (strcmp (vodeoptions.NormControl, 'on')) vnormcontrol = true;
else vnormcontrol = false; end
%# Implementation of the option NonNegative has been finished. This
%# option can be set by the user to another value than default value.
if (~isempty (vodeoptions.NonNegative))
if (isempty (vodeoptions.Mass)), vhavenonnegative = true;
else
vhavenonnegative = false;
warning ('OdePkg:InvalidArgument', ...
'Option "NonNegative" will be ignored if mass matrix is set');
end
else vhavenonnegative = false;
end
%# Implementation of the option OutputFcn has been finished. This
%# option can be set by the user to another value than default value.
if (isempty (vodeoptions.OutputFcn) && nargout == 0)
vodeoptions.OutputFcn = @odeplot;
vhaveoutputfunction = true;
elseif (isempty (vodeoptions.OutputFcn)), vhaveoutputfunction = false;
else vhaveoutputfunction = true;
end
%# Implementation of the option OutputSel has been finished. This
%# option can be set by the user to another value than default value.
if (~isempty (vodeoptions.OutputSel)), vhaveoutputselection = true;
else vhaveoutputselection = false; end
%# Implementation of the option OutputSave has been finished. This
%# option can be set by the user to another value than default value.
if (isempty (vodeoptions.OutputSave)), vodeoptions.OutputSave = 1;
end
%# Implementation of the option Refine has been finished. This option
%# can be set by the user to another value than default value.
if (vodeoptions.Refine > 0), vhaverefine = true;
else vhaverefine = false; end
%# Implementation of the option Stats has been finished. This option
%# can be set by the user to another value than default value.
%# Implementation of the option InitialStep has been finished. This
%# option can be set by the user to another value than default value.
if (isempty (vodeoptions.InitialStep) && ~vstepsizefixed)
vodeoptions.InitialStep = (vslot(1,2) - vslot(1,1)) / 10;
vodeoptions.InitialStep = vodeoptions.InitialStep / 10^vodeoptions.Refine;
warning ('OdePkg:InvalidArgument', ...
'Option "InitialStep" not set, new value %f is used', vodeoptions.InitialStep);
end
%# Implementation of the option MaxStep has been finished. This option
%# can be set by the user to another value than default value.
if (isempty (vodeoptions.MaxStep) && ~vstepsizefixed)
vodeoptions.MaxStep = abs (vslot(1,2) - vslot(1,1)) / 10;
warning ('OdePkg:InvalidArgument', ...
'Option "MaxStep" not set, new value %f is used', vodeoptions.MaxStep);
end
%# Implementation of the option Events has been finished. This option
%# can be set by the user to another value than default value.
if (~isempty (vodeoptions.Events)), vhaveeventfunction = true;
else vhaveeventfunction = false; end
%# The options 'Jacobian', 'JPattern' and 'Vectorized' will be ignored
%# by this solver because this solver uses an explicit Runge-Kutta
%# method and therefore no Jacobian calculation is necessary
if (~isequal (vodeoptions.Jacobian, vodetemp.Jacobian))
warning ('OdePkg:InvalidArgument', ...
'Option "Jacobian" will be ignored by this solver');
end
if (~isequal (vodeoptions.JPattern, vodetemp.JPattern))
warning ('OdePkg:InvalidArgument', ...
'Option "JPattern" will be ignored by this solver');
end
if (~isequal (vodeoptions.Vectorized, vodetemp.Vectorized))
warning ('OdePkg:InvalidArgument', ...
'Option "Vectorized" will be ignored by this solver');
end
if (~isequal (vodeoptions.NewtonTol, vodetemp.NewtonTol))
warning ('OdePkg:InvalidArgument', ...
'Option "NewtonTol" will be ignored by this solver');
end
if (~isequal (vodeoptions.MaxNewtonIterations,...
vodetemp.MaxNewtonIterations))
warning ('OdePkg:InvalidArgument', ...
'Option "MaxNewtonIterations" will be ignored by this solver');
end
%# Implementation of the option Mass has been finished. This option
%# can be set by the user to another value than default value.
if (~isempty (vodeoptions.Mass) && isnumeric (vodeoptions.Mass))
vhavemasshandle = false; vmass = vodeoptions.Mass; %# constant mass
elseif (isa (vodeoptions.Mass, 'function_handle'))
vhavemasshandle = true; %# mass defined by a function handle
else %# no mass matrix - creating a diag-matrix of ones for mass
vhavemasshandle = false; %# vmass = diag (ones (length (vinit), 1), 0);
end
%# Implementation of the option MStateDependence has been finished.
%# This option can be set by the user to another value than default
%# value.
if (strcmp (vodeoptions.MStateDependence, 'none'))
vmassdependence = false;
else vmassdependence = true;
end
%# Other options that are not used by this solver. Print a warning
%# message to tell the user that the option(s) is/are ignored.
if (~isequal (vodeoptions.MvPattern, vodetemp.MvPattern))
warning ('OdePkg:InvalidArgument', ...
'Option "MvPattern" will be ignored by this solver');
end
if (~isequal (vodeoptions.MassSingular, vodetemp.MassSingular))
warning ('OdePkg:InvalidArgument', ...
'Option "MassSingular" will be ignored by this solver');
end
if (~isequal (vodeoptions.InitialSlope, vodetemp.InitialSlope))
warning ('OdePkg:InvalidArgument', ...
'Option "InitialSlope" will be ignored by this solver');
end
if (~isequal (vodeoptions.MaxOrder, vodetemp.MaxOrder))
warning ('OdePkg:InvalidArgument', ...
'Option "MaxOrder" will be ignored by this solver');
end
if (~isequal (vodeoptions.BDF, vodetemp.BDF))
warning ('OdePkg:InvalidArgument', ...
'Option "BDF" will be ignored by this solver');
end
%# Starting the initialisation of the core solver ode23
vtimestamp = vslot(1,1); %# timestamp = start time
vtimelength = length (vslot); %# length needed if fixed steps
vtimestop = vslot(1,vtimelength); %# stop time = last value
%# 20110611, reported by Nils Strunk
%# Make it possible to solve equations from negativ to zero,
%# eg. vres = ode23 (@(t,y) y, [-2 0], 2);
vdirection = sign (vtimestop - vtimestamp); %# Direction flag
if (~vstepsizefixed)
if (sign (vodeoptions.InitialStep) == vdirection)
vstepsize = vodeoptions.InitialStep;
else %# Fix wrong direction of InitialStep.
vstepsize = - vodeoptions.InitialStep;
end
vminstepsize = (vtimestop - vtimestamp) / (1/eps);
else %# If step size is given then use the fixed time steps
vstepsize = vslot(1,2) - vslot(1,1);
vminstepsize = sign (vstepsize) * eps;
end
vretvaltime = vtimestamp; %# first timestamp output
vretvalresult = vinit; %# first solution output
%# Initialize the OutputFcn
if (vhaveoutputfunction)
if (vhaveoutputselection) vretout = vretvalresult(vodeoptions.OutputSel);
else vretout = vretvalresult; end
feval (vodeoptions.OutputFcn, vslot.', ...
vretout.', 'init', vfunarguments{:});
end
%# Initialize the EventFcn
if (vhaveeventfunction)
odepkg_event_handle (vodeoptions.Events, vtimestamp, ...
vretvalresult.', 'init', vfunarguments{:});
end
vpow = 1/3; %# 20071016, reported by Luis Randez
va = [ 0, 0, 0; %# The Runge-Kutta-Fehlberg 2(3) coefficients
1/2, 0, 0; %# Coefficients proved on 20060827
-1, 2, 0]; %# See p.91 in Ascher & Petzold
vb2 = [0; 1; 0]; %# 2nd and 3rd order
vb3 = [1/6; 2/3; 1/6]; %# b-coefficients
vc = sum (va, 2);
%# The solver main loop - stop if the endpoint has been reached
vcntloop = 2; vcntcycles = 1; vu = vinit; vk = vu.' * zeros(1,3);
vcntiter = 0; vunhandledtermination = true; vcntsave = 2;
while ((vdirection * (vtimestamp) < vdirection * (vtimestop)) && ...
(vdirection * (vstepsize) >= vdirection * (vminstepsize)))
%# Hit the endpoint of the time slot exactely
if (vdirection * (vtimestamp + vstepsize) > vdirection * vtimestop)
%# vstepsize = vtimestop - vdirection * vtimestamp;
%# 20110611, reported by Nils Strunk
%# The endpoint of the time slot must be hit exactly,
%# eg. vsol = ode23 (@(t,y) y, [0 -1], 1);
vstepsize = vdirection * abs (abs (vtimestop) - abs (vtimestamp));
end
%# Estimate the three results when using this solver
for j = 1:3
vthetime = vtimestamp + vc(j,1) * vstepsize;
vtheinput = vu.' + vstepsize * vk(:,1:j-1) * va(j,1:j-1).';
if (vhavemasshandle) %# Handle only the dynamic mass matrix,
if (vmassdependence) %# constant mass matrices have already
vmass = feval ... %# been set before (if any)
(vodeoptions.Mass, vthetime, vtheinput, vfunarguments{:});
else %# if (vmassdependence == false)
vmass = feval ... %# then we only have the time argument
(vodeoptions.Mass, vthetime, vfunarguments{:});
end
vk(:,j) = vmass \ feval ...
(vfun, vthetime, vtheinput, vfunarguments{:});
else
vk(:,j) = feval ...
(vfun, vthetime, vtheinput, vfunarguments{:});
end
end
%# Compute the 2nd and the 3rd order estimation
y2 = vu.' + vstepsize * (vk * vb2);
y3 = vu.' + vstepsize * (vk * vb3);
if (vhavenonnegative)
vu(vodeoptions.NonNegative) = abs (vu(vodeoptions.NonNegative));
y2(vodeoptions.NonNegative) = abs (y2(vodeoptions.NonNegative));
y3(vodeoptions.NonNegative) = abs (y3(vodeoptions.NonNegative));
end
if (vhaveoutputfunction && vhaverefine)
vSaveVUForRefine = vu;
end
%# Calculate the absolute local truncation error and the acceptable error
if (~vstepsizefixed)
if (~vnormcontrol)
vdelta = abs (y3 - y2);
vtau = max (vodeoptions.RelTol * abs (vu.'), vodeoptions.AbsTol);
else
vdelta = norm (y3 - y2, Inf);
vtau = max (vodeoptions.RelTol * max (norm (vu.', Inf), 1.0), ...
vodeoptions.AbsTol);
end
else %# if (vstepsizefixed == true)
vdelta = 1; vtau = 2;
end
%# If the error is acceptable then update the vretval variables
if (all (vdelta <= vtau))
vtimestamp = vtimestamp + vstepsize;
vu = y3.'; %# MC2001: the higher order estimation as "local extrapolation"
%# Save the solution every vodeoptions.OutputSave steps
if (mod (vcntloop-1,vodeoptions.OutputSave) == 0)
vretvaltime(vcntsave,:) = vtimestamp;
vretvalresult(vcntsave,:) = vu;
vcntsave = vcntsave + 1;
end
vcntloop = vcntloop + 1; vcntiter = 0;
%# Call plot only if a valid result has been found, therefore this
%# code fragment has moved here. Stop integration if plot function
%# returns false
if (vhaveoutputfunction)
for vcnt = 0:vodeoptions.Refine %# Approximation between told and t
if (vhaverefine) %# Do interpolation
vapproxtime = (vcnt + 1) * vstepsize / (vodeoptions.Refine + 2);
vapproxvals = vSaveVUForRefine.' + vapproxtime * (vk * vb3);
vapproxtime = (vtimestamp - vstepsize) + vapproxtime;
else
vapproxvals = vu.';
vapproxtime = vtimestamp;
end
if (vhaveoutputselection)
vapproxvals = vapproxvals(vodeoptions.OutputSel);
end
vpltret = feval (vodeoptions.OutputFcn, vapproxtime, ...
vapproxvals, [], vfunarguments{:});
if vpltret %# Leave refinement loop
break;
end
end
if (vpltret) %# Leave main loop
vunhandledtermination = false;
break;
end
end
%# Call event only if a valid result has been found, therefore this
%# code fragment has moved here. Stop integration if veventbreak is
%# true
if (vhaveeventfunction)
vevent = ...
odepkg_event_handle (vodeoptions.Events, vtimestamp, ...
vu(:), [], vfunarguments{:});
if (~isempty (vevent{1}) && vevent{1} == 1)
vretvaltime(vcntloop-1,:) = vevent{3}(end,:);
vretvalresult(vcntloop-1,:) = vevent{4}(end,:);
vunhandledtermination = false; break;
end
end
end %# If the error is acceptable ...
%# Update the step size for the next integration step
if (~vstepsizefixed)
%# 20080425, reported by Marco Caliari
%# vdelta cannot be negative (because of the absolute value that
%# has been introduced) but it could be 0, then replace the zeros
%# with the maximum value of vdelta
vdelta(find (vdelta == 0)) = max (vdelta);
%# It could happen that max (vdelta) == 0 (ie. that the original
%# vdelta was 0), in that case we double the previous vstepsize
vdelta(find (vdelta == 0)) = max (vtau) .* (0.4 ^ (1 / vpow));
if (vdirection == 1)
vstepsize = min (vodeoptions.MaxStep, ...
min (0.8 * vstepsize * (vtau ./ vdelta) .^ vpow));
else
vstepsize = max (- vodeoptions.MaxStep, ...
max (0.8 * vstepsize * (vtau ./ vdelta) .^ vpow));
end
else %# if (vstepsizefixed)
if (vcntloop <= vtimelength)
vstepsize = vslot(vcntloop) - vslot(vcntloop-1);
else %# Get out of the main integration loop
break;
end
end
%# Update counters that count the number of iteration cycles
vcntcycles = vcntcycles + 1; %# Needed for cost statistics
vcntiter = vcntiter + 1; %# Needed to find iteration problems
%# Stop solving because the last 1000 steps no successful valid
%# value has been found
if (vcntiter >= 5000)
error (['Solving has not been successful. The iterative', ...
' integration loop exited at time t = %f before endpoint at', ...
' tend = %f was reached. This happened because the iterative', ...
' integration loop does not find a valid solution at this time', ...
' stamp. Try to reduce the value of "InitialStep" and/or', ...
' "MaxStep" with the command "odeset".\n'], vtimestamp, vtimestop);
end
end %# The main loop
%# Check if integration of the ode has been successful
if (vdirection * vtimestamp < vdirection * vtimestop)
if (vunhandledtermination == true)
error ('OdePkg:InvalidArgument', ...
['Solving has not been successful. The iterative', ...
' integration loop exited at time t = %f', ...
' before endpoint at tend = %f was reached. This may', ...
' happen if the stepsize grows smaller than defined in', ...
' vminstepsize. Try to reduce the value of "InitialStep" and/or', ...
' "MaxStep" with the command "odeset".\n'], vtimestamp, vtimestop);
else
warning ('OdePkg:InvalidArgument', ...
['Solver has been stopped by a call of "break" in', ...
' the main iteration loop at time t = %f before endpoint at', ...
' tend = %f was reached. This may happen because the @odeplot', ...
' function returned "true" or the @event function returned "true".'], ...
vtimestamp, vtimestop);
end
end
%# Postprocessing, do whatever when terminating integration algorithm
if (vhaveoutputfunction) %# Cleanup plotter
feval (vodeoptions.OutputFcn, vtimestamp, ...
vu.', 'done', vfunarguments{:});
end
if (vhaveeventfunction) %# Cleanup event function handling
odepkg_event_handle (vodeoptions.Events, vtimestamp, ...
vu.', 'done', vfunarguments{:});
end
%# Save the last step, if not already saved
if (mod (vcntloop-2,vodeoptions.OutputSave) ~= 0)
vretvaltime(vcntsave,:) = vtimestamp;
vretvalresult(vcntsave,:) = vu;
end
%# Print additional information if option Stats is set
if (strcmp (vodeoptions.Stats, 'on'))
vhavestats = true;
vnsteps = vcntloop-2; %# vcntloop from 2..end
vnfailed = (vcntcycles-1)-(vcntloop-2)+1; %# vcntcycl from 1..end
vnfevals = 3*(vcntcycles-1); %# number of ode evaluations
vndecomps = 0; %# number of LU decompositions
vnpds = 0; %# number of partial derivatives
vnlinsols = 0; %# no. of solutions of linear systems
%# Print cost statistics if no output argument is given
if (nargout == 0)
vmsg = fprintf (1, 'Number of successful steps: %d\n', vnsteps);
vmsg = fprintf (1, 'Number of failed attempts: %d\n', vnfailed);
vmsg = fprintf (1, 'Number of function calls: %d\n', vnfevals);
end
else
vhavestats = false;
end
if (nargout == 1) %# Sort output variables, depends on nargout
varargout{1}.x = vretvaltime; %# Time stamps are saved in field x
varargout{1}.y = vretvalresult; %# Results are saved in field y
varargout{1}.solver = 'ode23'; %# Solver name is saved in field solver
if (vhaveeventfunction)
varargout{1}.ie = vevent{2}; %# Index info which event occured
varargout{1}.xe = vevent{3}; %# Time info when an event occured
varargout{1}.ye = vevent{4}; %# Results when an event occured
end
if (vhavestats)
varargout{1}.stats = struct;
varargout{1}.stats.nsteps = vnsteps;
varargout{1}.stats.nfailed = vnfailed;
varargout{1}.stats.nfevals = vnfevals;
varargout{1}.stats.npds = vnpds;
varargout{1}.stats.ndecomps = vndecomps;
varargout{1}.stats.nlinsols = vnlinsols;
end
elseif (nargout == 2)
varargout{1} = vretvaltime; %# Time stamps are first output argument
varargout{2} = vretvalresult; %# Results are second output argument
elseif (nargout == 5)
varargout{1} = vretvaltime; %# Same as (nargout == 2)
varargout{2} = vretvalresult; %# Same as (nargout == 2)
varargout{3} = []; %# LabMat doesn't accept lines like
varargout{4} = []; %# varargout{3} = varargout{4} = [];
varargout{5} = [];
if (vhaveeventfunction)
varargout{3} = vevent{3}; %# Time info when an event occured
varargout{4} = vevent{4}; %# Results when an event occured
varargout{5} = vevent{2}; %# Index info which event occured
end
end
end
%! # We are using the "Van der Pol" implementation for all tests that
%! # are done for this function. We also define a Jacobian, Events,
%! # pseudo-Mass implementation. For further tests we also define a
%! # reference solution (computed at high accuracy) and an OutputFcn
%!function [ydot] = fpol (vt, vy, varargin) %# The Van der Pol
%! ydot = [vy(2); (1 - vy(1)^2) * vy(2) - vy(1)];
%!function [vjac] = fjac (vt, vy, varargin) %# its Jacobian
%! vjac = [0, 1; -1 - 2 * vy(1) * vy(2), 1 - vy(1)^2];
%!function [vjac] = fjcc (vt, vy, varargin) %# sparse type
%! vjac = sparse ([0, 1; -1 - 2 * vy(1) * vy(2), 1 - vy(1)^2]);
%!function [vval, vtrm, vdir] = feve (vt, vy, varargin)
%! vval = fpol (vt, vy, varargin); %# We use the derivatives
%! vtrm = zeros (2,1); %# that's why component 2
%! vdir = ones (2,1); %# seems to not be exact
%!function [vval, vtrm, vdir] = fevn (vt, vy, varargin)
%! vval = fpol (vt, vy, varargin); %# We use the derivatives
%! vtrm = ones (2,1); %# that's why component 2
%! vdir = ones (2,1); %# seems to not be exact
%!function [vmas] = fmas (vt, vy)
%! vmas = [1, 0; 0, 1]; %# Dummy mass matrix for tests
%!function [vmas] = fmsa (vt, vy)
%! vmas = sparse ([1, 0; 0, 1]); %# A sparse dummy matrix
%!function [vref] = fref () %# The computed reference sol
%! vref = [0.32331666704577, -1.83297456798624];
%!function [vout] = fout (vt, vy, vflag, varargin)
%! if (regexp (char (vflag), 'init') == 1)
%! if (any (size (vt) ~= [2, 1])) error ('"fout" step "init"'); end
%! elseif (isempty (vflag))
%! if (any (size (vt) ~= [1, 1])) error ('"fout" step "calc"'); end
%! vout = false;
%! elseif (regexp (char (vflag), 'done') == 1)
%! if (any (size (vt) ~= [1, 1])) error ('"fout" step "done"'); end
%! else error ('"fout" invalid vflag');
%! end
%!
%! %# Turn off output of warning messages for all tests, turn them on
%! %# again if the last test is called
%!error %# input argument number one
%! warning ('off', 'OdePkg:InvalidArgument');
%! B = ode23 (1, [0 25], [3 15 1]);
%!error %# input argument number two
%! B = ode23 (@fpol, 1, [3 15 1]);
%!error %# input argument number three
%! B = ode23 (@flor, [0 25], 1);
%!test %# one output argument
%! vsol = ode23 (@fpol, [0 2], [2 0]);
%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
%! assert (isfield (vsol, 'solver'));
%! assert (vsol.solver, 'ode23');
%!test %# two output arguments
%! [vt, vy] = ode23 (@fpol, [0 2], [2 0]);
%! assert ([vt(end), vy(end,:)], [2, fref], 1e-3);
%!test %# five output arguments and no Events
%! [vt, vy, vxe, vye, vie] = ode23 (@fpol, [0 2], [2 0]);
%! assert ([vt(end), vy(end,:)], [2, fref], 1e-3);
%! assert ([vie, vxe, vye], []);
%!test %# anonymous function instead of real function
%! fvdb = @(vt,vy) [vy(2); (1 - vy(1)^2) * vy(2) - vy(1)];
%! vsol = ode23 (fvdb, [0 2], [2 0]);
%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
%!test %# extra input arguments passed through
%! vsol = ode23 (@fpol, [0 2], [2 0], 12, 13, 'KL');
%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
%!test %# empty OdePkg structure *but* extra input arguments
%! vopt = odeset;
%! vsol = ode23 (@fpol, [0 2], [2 0], vopt, 12, 13, 'KL');
%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
%!error %# strange OdePkg structure
%! vopt = struct ('foo', 1);
%! vsol = ode23 (@fpol, [0 2], [2 0], vopt);
%!test %# Solve vdp in fixed step sizes
%! vsol = ode23 (@fpol, [0:0.1:2], [2 0]);
%! assert (vsol.x(:), [0:0.1:2]');
%! assert (vsol.y(end,:), fref, 1e-3);
%!test %# Solve in backward direction starting at t=0
%! vref = [-1.205364552835178, 0.951542399860817];
%! vsol = ode23 (@fpol, [0 -2], [2 0]);
%! assert ([vsol.x(end), vsol.y(end,:)], [-2, vref], 1e-3);
%!test %# Solve in backward direction starting at t=2
%! vref = [-1.205364552835178, 0.951542399860817];
%! vsol = ode23 (@fpol, [2 -2], fref);
%! assert ([vsol.x(end), vsol.y(end,:)], [-2, vref], 1e-3);
%!test %# Solve another anonymous function in backward direction
%! vref = [-1, 0.367879437558975];
%! vsol = ode23 (@(t,y) y, [0 -1], 1);
%! assert ([vsol.x(end), vsol.y(end,:)], vref, 1e-3);
%!test %# Solve another anonymous function below zero
%! vref = [0, 14.77810590694212];
%! vsol = ode23 (@(t,y) y, [-2 0], 2);
%! assert ([vsol.x(end), vsol.y(end,:)], vref, 1e-3);
%!test %# AbsTol option
%! vopt = odeset ('AbsTol', 1e-5);
%! vsol = ode23 (@fpol, [0 2], [2 0], vopt);
%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
%!test %# AbsTol and RelTol option
%! vopt = odeset ('AbsTol', 1e-8, 'RelTol', 1e-8);
%! vsol = ode23 (@fpol, [0 2], [2 0], vopt);
%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
%!test %# RelTol and NormControl option -- higher accuracy
%! vopt = odeset ('RelTol', 1e-8, 'NormControl', 'on');
%! vsol = ode23 (@fpol, [0 2], [2 0], vopt);
%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-6);
%!test %# Keeps initial values while integrating
%! vopt = odeset ('NonNegative', 2);
%! vsol = ode23 (@fpol, [0 2], [2 0], vopt);
%! assert ([vsol.x(end), vsol.y(end,:)], [2, 2, 0], 1e-1);
%!test %# Details of OutputSel and Refine can't be tested
%! vopt = odeset ('OutputFcn', @fout, 'OutputSel', 1, 'Refine', 5);
%! vsol = ode23 (@fpol, [0 2], [2 0], vopt);
%!test %# Details of OutputSave can't be tested
%! vopt = odeset ('OutputSave', 1, 'OutputSel', 1);
%! vsla = ode23 (@fpol, [0 2], [2 0], vopt);
%! vopt = odeset ('OutputSave', 2);
%! vslb = ode23 (@fpol, [0 2], [2 0], vopt);
%! assert (length (vsla.x) > length (vslb.x))
%!test %# Stats must add further elements in vsol
%! vopt = odeset ('Stats', 'on');
%! vsol = ode23 (@fpol, [0 2], [2 0], vopt);
%! assert (isfield (vsol, 'stats'));
%! assert (isfield (vsol.stats, 'nsteps'));
%!test %# InitialStep option
%! vopt = odeset ('InitialStep', 1e-8);
%! vsol = ode23 (@fpol, [0 0.2], [2 0], vopt);
%! assert ([vsol.x(2)-vsol.x(1)], [1e-8], 1e-9);
%!test %# MaxStep option
%! vopt = odeset ('MaxStep', 1e-2);
%! vsol = ode23 (@fpol, [0 0.2], [2 0], vopt);
%! assert ([vsol.x(5)-vsol.x(4)], [1e-2], 1e-2);
%!test %# Events option add further elements in vsol
%! vopt = odeset ('Events', @feve);
%! vsol = ode23 (@fpol, [0 10], [2 0], vopt);
%! assert (isfield (vsol, 'ie'));
%! assert (vsol.ie, [2; 1; 2; 1]);
%! assert (isfield (vsol, 'xe'));
%! assert (isfield (vsol, 'ye'));
%!test %# Events option, now stop integration
%! vopt = odeset ('Events', @fevn, 'NormControl', 'on');
%! vsol = ode23 (@fpol, [0 10], [2 0], vopt);
%! assert ([vsol.ie, vsol.xe, vsol.ye], ...
%! [2.0, 2.496110, -0.830550, -2.677589], 1e-3);
%!test %# Events option, five output arguments
%! vopt = odeset ('Events', @fevn, 'NormControl', 'on');
%! [vt, vy, vxe, vye, vie] = ode23 (@fpol, [0 10], [2 0], vopt);
%! assert ([vie, vxe, vye], ...
%! [2.0, 2.496110, -0.830550, -2.677589], 1e-3);
%!test %# Jacobian option
%! vopt = odeset ('Jacobian', @fjac);
%! vsol = ode23 (@fpol, [0 2], [2 0], vopt);
%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
%!test %# Jacobian option and sparse return value
%! vopt = odeset ('Jacobian', @fjcc);
%! vsol = ode23 (@fpol, [0 2], [2 0], vopt);
%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
%!
%! %# test for JPattern option is missing
%! %# test for Vectorized option is missing
%! %# test for NewtonTol option is missing
%! %# test for MaxNewtonIterations option is missing
%!
%!test %# Mass option as function
%! vopt = odeset ('Mass', @fmas);
%! vsol = ode23 (@fpol, [0 2], [2 0], vopt);
%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
%!test %# Mass option as matrix
%! vopt = odeset ('Mass', eye (2,2));
%! vsol = ode23 (@fpol, [0 2], [2 0], vopt);
%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
%!test %# Mass option as sparse matrix
%! vopt = odeset ('Mass', sparse (eye (2,2)));
%! vsol = ode23 (@fpol, [0 2], [2 0], vopt);
%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
%!test %# Mass option as function and sparse matrix
%! vopt = odeset ('Mass', @fmsa);
%! vsol = ode23 (@fpol, [0 2], [2 0], vopt);
%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
%!test %# Mass option as function and MStateDependence
%! vopt = odeset ('Mass', @fmas, 'MStateDependence', 'strong');
%! vsol = ode23 (@fpol, [0 2], [2 0], vopt);
%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
%!test %# Set BDF option to something else than default
%! vopt = odeset ('BDF', 'on');
%! [vt, vy] = ode23 (@fpol, [0 2], [2 0], vopt);
%! assert ([vt(end), vy(end,:)], [2, fref], 1e-3);
%!
%! %# test for MvPattern option is missing
%! %# test for InitialSlope option is missing
%! %# test for MaxOrder option is missing
%!
%! warning ('on', 'OdePkg:InvalidArgument');
%# Local Variables: ***
%# mode: octave ***
%# End: ***
\showoff
|
Profil
|
Delastelle
Senior
Dabei seit: 17.11.2006
Mitteilungen: 2320
| Beitrag No.2, eingetragen 2020-03-04
|
Hallo Specialagent!
Du hast das Thema abgehakt.
Mich würde nur interessieren ob der Tipp mit der veränderten Toleranz (atol) zum Ziel geführt hat!
Wenn eine numerische Zahl mit Typ double eine Genaugigkeit von ca.15 bis 16 Stellen hat, kann man normalerweise keine Genauigkeit von 1e-13 in einem Verfahren erreichen.
Viele Grüße
Ronald
|
Profil
|
Specialagent hat die Antworten auf ihre/seine Frage gesehen.
Specialagent hatte hier bereits selbst das Ok-Häkchen gesetzt. |
|
All logos and trademarks in this site are property of their respective owner. The comments are property of their posters, all the rest © 2001-2023 by Matroids Matheplanet
This web site was originally made with PHP-Nuke, a former web portal system written in PHP that seems no longer to be maintained nor supported. PHP-Nuke is Free Software released under the GNU/GPL license.
Ich distanziere mich von rechtswidrigen oder anstößigen Inhalten, die sich trotz aufmerksamer Prüfung hinter hier verwendeten Links verbergen mögen.
Lesen Sie die
Nutzungsbedingungen,
die Distanzierung,
die Datenschutzerklärung und das Impressum.
[Seitenanfang]
|
Home / Seite ohne Frame
Wie unterstütze ich den Matheplaneten mit Geld?
