function [RESULT, SETTINGS]=ivp_solve(SETTINGS)
% function thal solves the initial value problem
% by using the specified parameters out of SETTINGS
%
% input: o SETTINGS....struct containing all information about the
% current program run
%
% example o SETTINGS =
%
% step: 0.1000
% method: {'explicit Euler-method'}
% calculate_region_of_stability: 0
% equation: @(x,y)(cos(y)).^2./sin(y).*x.*exp(1/2.*x)
% fu_exact: @(x)acos(1/(4*exp(1/2*x)*(1/2*x-1)+1/cos(pi/2-0.2626)*(14*exp(-5/2)*cos(pi/2-0.2626)+1)))
% beginvalue: -5
% endvalue: 4
% startvalue: 1.3082
% isRKV: 1
% step_control: 0
%
% output: o RESULT....struct containing all needed information about the
% numerical soltuion
%
% example: o RESULT =
%
% numerical_solution: [1x91 double]
% exact_solution: [1x91 double]
% gridpoints: [1x91 double]
%
% notes: o
%
% authors: o Christian Jkel (University of Technology Dresden)
% o Daniel Klawitter (University of Technology Dresden)
stop_calculation = 0;
f = SETTINGS.equation;
% calculates the k's for implicit RKV
function [ka] = k_imp(j,age)
function[Funktion] = F(x,y1,y2)
Funktion = y2-f(x+a*age,(y1'+age*B*y2')');
end
ka = fsolve(@(y)F(step_vec(j),u(j),y),ones(1,s),optimset('Display','off'));
end
% calculates the k's for explicit RKV
function [ka] = k(j,age)
for i2 = 1:s
ka(1) = f(step_vec(j)+a(i2)*age,u(j));
ka(i2) = f(step_vec(j)+a(i2)*age,u(j)+age*sum(B(i2,1:i2-1).*ka(1:i2-1)));
end
end
% this function calculates initial values for the linear multistep methods
function[startwerte] = berechne_Startwerte(stufe)
a = [0, 1/5, 3/10, 4/5, 8/9, 1, 1];
B = [0 , 0, 0, 0, 0, 0, 0;
1/5 , 0, 0, 0, 0, 0, 0;
3/40 , 9/40, 0, 0, 0, 0, 0;
44/45 , -56/15, 32/9, 0, 0, 0, 0;
19372/6561, -25360/2187, 64448/6561, -212/729, 0, 0, 0;
9017/3168 , -355/33, 46732/5247, 49/176, -5103/18656, 0, 0;
35/384 , 0, 500/1113, 125/192, -2187/6784, 11/84, 0];
c = [5179/57600,0,7571/16695,393/640,-92097/339200,187/2100,1/40];
s = 7;
p=5;
tele(1)=SETTINGS.startvalue;
function [ka] = kaa(j,age)
for i2 = 1:s
ka(1) = f(step_vec(j)+a(i2)*age,tele(j));
ka(i2) = f(step_vec(j)+a(i2)*age,tele(j)+age*sum(B(i2,1:i2-1).*ka(1:i2-1)));
end
end
if SETTINGS.step_control==1
i=2;
method=@(sta,z,vector)(vector(z)+sta*sum(c.*kaa(z,sta)));
while go==true
u1 = method(h_old,i-1,tele);
u2 = method(h_old/2,i-1,tele); % extrapolation to estimate the error
est=abs((2^p*u2-u1)/(2^p-1)-u1);
if est<TOL & est>0.1*TOL % step size will be accepted
RESULT.steps(i)=h_old;
step_vec(i)=step_vec(i-1)+h_old;
tele(i)=u1;
i=i+1;
else
h_new=h_old*(0.5*TOL/est)^(1/(p+1));
if h_new<h_old % decrease of the step size
h_old=h_new;
elseif h_new>=h_old & est<0.1*TOL % increase of the step size for the next step
tele(i)=u1;
step_vec(i)=step_vec(i-1)+h_old;
RESULT.steps(i)=h_old;
h_old=h_new;
i=i+1;
end
end
if i==length(rho)
go = false;
end
end
else
for i = 2:stufe
tele(i) = tele(i-1)+SETTINGS.step*sum(c.*kaa(i-1,SETTINGS.step));
end
end
startwerte = tele;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%% THE METHODS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch char(SETTINGS.method)
case {'explicit Euler-method'}
% Euler explizit --> explizit
B = [0];
a = [0];
c = [1];
s = 1;
p=1;
case {'Heun-Verfahren'}
% Heun-Verfahren --> explizit
B = [0, 0;
1, 0];
a = [0 , 1 ];
c = [1/2, 1/2];
s = 2;
p=2;
case {'Klassisches RKV'}
% klassisches RKV --> explizit
B = [0 , 0 , 0, 0;
1/2, 0 , 0, 0;
0 ,1/2, 0, 0;
0 , 0, 1, 0];
a = [0 , 1/2, 1/2 , 1];
c = [1/6, 1/3, 1/3, 1/6];
s = 4;
p=4;
case {'Dormand-Prince RKV'}
% embedded RKV from Dormand-Prince sheme
a = [0, 1/5, 3/10 , 4/5, 8/9, 1 , 1 ];
B = [0 , 0 , 0 , 0 , 0 , 0 , 0;
1/5 , 0 , 0 , 0 , 0 , 0 , 0;
3/40 , 9/40 , 0 , 0 , 0 , 0 , 0;
44/45 , -56/15 , 32/9 , 0 , 0 , 0 , 0;
19372/6561, -25360/2187, 64448/6561 , -212/729, 0 , 0 , 0;
9017/3168 , -355/33 , 46732/5247 , 49/176 , -5103/18656 , 0 , 0;
35/384 , 0 , 500/1113 , 125/192 , -2187/6784 , 11/84 , 0];
% with this c order4
c = [5179/57600 , 0, 7571/16695 , 393/640, -92097/339200, 187/2100, 1/40];
% and with this one order 5
c = [35/384 , 0 , 500/1113 , 125/192 , -2187/6784 , 11/84 , 0];
s = 7;
p=5;
case {'implicit Euler-method'}
% implizites verfahren
B = [1];
a = [1];
c = [1];
s = 1;
p=1;
case {'Radau II RKV'}
% Radau IIA method order 3 entnommen von
% http://www-user.tu-chemnitz.de/~benner/Lehre/NumerikODE/Folie_IRK.pdf
B = [5/12, -1/12;
3/4 , 1/4 ];
a = [1/3,1];
c = [3/4,1/4];
s = 2;
p=3;
case {'Gau RKV'}
% gauverfahren der ordnung 6 entnommen von
% http://www-user.tu-chemnitz.de/~benner/Lehre/NumerikODE/Folie_IRK.pdf
B = [5/36 , (10-3*sqrt(15))/45, (25-6*sqrt(15))/180;
(10+3*sqrt(15))/72 , 2/9 , (10-3*sqrt(15))/72;
(25+6*sqrt(15))/180, (10+3*sqrt(15))/45, 5/36];
a = [(5-sqrt(15))/10,1/2,(5+sqrt(15))/10];
c = [5/18,4/9,5/18];
s = 3;
p=6;
end
%%%%%%%%%%%%%%%%% calculation of the numerical solution %%%%%%%%%%%%%%%%%%%%%%
step_vec(1)=SETTINGS.beginvalue; % vector with grid points
h_old=SETTINGS.step;
if SETTINGS.step_control==1
RESULT.steps(1)=SETTINGS.step; % vektor with step size for each grid point
go=true;
TOL=SETTINGS.tolerance;
else
step_vec = SETTINGS.beginvalue:SETTINGS.step:SETTINGS.endvalue;
end
switch char(SETTINGS.method)
case {'explicit Euler-method','Heun-Verfahren','Klassisches RKV','Dormand-Prince RKV'}
% explizite RKV
u(1) = SETTINGS.startvalue;
method=@(sta,z,vector,vector_step)(vector(z)+sta*sum(c.*k(z,sta)));
i=2;
case {'implicit Euler-method','Radau II RKV','Gau RKV','radauIIA'}
u(1) = SETTINGS.startvalue;
method=@(sta,z,vector,vector_step)(vector(z)+sta*sum(c.*k_imp(z,sta)));
i=2;
case {'Adams-Bashforth(2)'}
% ein zweistufiges Adamsverfahren (explizit)
rho = [1, -1 , 0 ];
sigma = [0, 3/2, -1/2];
u(1) = SETTINGS.startvalue;
startwerte = berechne_Startwerte(2);
u(2) = startwerte(2);
method = @(sta,z,vector,vector_step)(vector(z)+sta/2*(3*f(vector_step(z),vector(z))-f(vector_step(z-1),vector(z-1))));
i=3;
go=true;
case {'Adams-Bashforth(4)'}
% ein vierstufiges Verfahren (explizit)
rho = [1, -1 , 0 , 0 ,0 ];
sigma = [0, 55/24, -59/24, 37/24,-9/24];
u(1) = SETTINGS.startvalue;
startwerte = berechne_Startwerte(4);
u(2:4) = startwerte(2:4);
i=5;
go=true;
method = @(sta,z,vector,vector_step)(vector(z)+sta/24*(55*f(vector_step(z),vector(z))...
-59*f(vector_step(z-1),vector(z-1))...
+37*f(vector_step(z-2),vector(z-2))...
-9*f(vector_step(z-3),vector(z-3))));
case {'Adams-Moulton(2)'}
% Adams-Moulton-verfahren, implizit, 2-stufig
rho = [1,-1,0];
sigma = [5/12,8/12,-1/12];
u(1) = SETTINGS.startvalue;
startwerte = berechne_Startwerte(2);
u(2) = startwerte(2);
i=3;
go=true;
method = @(sta,z,vector,vector_step)(fzero(@(y) y-vector(z)-sta/12*(5*f(vector_step(z+1),y)...
+8*f(vector_step(z),vector(z))...
-f(vector_step(z-1),vector(z-1))),vector(z)));
case {'Milne-Simpson(4)'}
% Milne-Simpson-verfahren, implizit, 4-stufig
rho = [1 , 0 , -1 , 0 , 0 ];
sigma = [29/90, 124/90, 24/90, 4/90, -1/90];
u(1) = SETTINGS.startvalue;
startwerte = berechne_Startwerte(4);
u(2:4) = startwerte(2:4);
i=5;
go=true;
method = @(sta,z,vector,vector_step) (fzero(@(y)y-vector(z-1)-sta/90*(29*f(vector_step(z+1),y)+...
124*f(vector_step(z),vector(z))...
+24*f(vector_step(z-1),vector(z-1))...
+4*f(vector_step(z-2),vector(z-2))...
-f(vector_step(z-3),vector(z-3))),vector(z)));
case {'BDF 5'}
% BDF-Verfahren, implizit, 5-stufig
rho = [137/60, -300/60, 300/60, -200/60, 75/60, -12/60];
sigma = [1 , 0 , 0 , 0 , 0 , 0];
u(1) = SETTINGS.startvalue;
startwerte = berechne_Startwerte(5);
u(2:5) = startwerte(2:5);
i=6;
go=true;
method = @(sta,z,vector,vector_step) (fzero(@(y) 1/60*(137*y-300*vector(z)+300*vector(z-1)-200*vector(z-2)+75*vector(z-3)-12*vector(z-4))...
-sta*f(vector_step(z+1),y),vector(z)));
case {'BDF 2'}
% BDF Verfahren 2-stufig, A-Stabil!!!!
rho=[3/2,-4/2,1/2];
sigma=[1,0,0];
u(1) = SETTINGS.startvalue;
startwerte = berechne_Startwerte(2);
u(2)=startwerte(2);
i=3;
go=true;
method = @(sta,z,vector,vector_step)(fzero(@(y)3/2*y-4/2*vector(z)+1/2*vector(z-1)-sta*f(vector_step(z+1),y),vector(z)));
case'eigenesRKV'
B = SETTINGS.B;
c = SETTINGS.c;
a = SETTINGS.a;
s = length(B);
u(1) = SETTINGS.startvalue;
i=2;
p=length(B); % rough estimation for the order of convergence.
method=@(sta,z,vector,vector_step)(vector(z)+sta*sum(c.*k_imp(z,sta)));
%if abs(u(i))>10^200
% errordlg('Your method produces a solution which can not be displayed. Change the chosen parameters!','error')
% stop_calculation=1;
% break
%
% end
case'eigenesMSV'
u(1) = SETTINGS.startvalue;
rho = SETTINGS.rho;
sigma = SETTINGS.sigma;
rho2 = fliplr(rho);
sigma2 = fliplr(sigma);
st = length(rho)-1;
startwerte = berechne_Startwerte(st);
u(2:st)= startwerte(2:st);
i=st+1;
go=true;
p=length(rho);
%zwischenwert1 = sum(rho2(1:st).*u(i-st:i-1));
%zwischenwert2 = sum(sigma2(1:st).*f(step_vec(i-st:i-1),u(i-st:i-1)));
method = @(sta,z,vector,vector_step)(fzero(@(y) y*rho(1)+sum(rho2(1:st).*vector(z+1-st:z))-sta*sigma(1)*f(vector_step(z+1),y)+sta*sum(sigma2(1:st).*f(vector_step(z+1-st:z),vector(z+1-st:z))),vector(z)));
%if abs(u(i))>10^200 | strcmp(num2str(u(i)),'NaN')
% errordlg('Your method produces a solution which can not be displayed. Change the chosen parameters!','error')
% stop_calculation=1;
% break
%
%end
end
if SETTINGS.step_control==1 % step control
while go==true
step_vec(i)=step_vec(i-1)+h_old; % this is for the implicit methods because they need this point for the error estimation
u1 = method(h_old,i-1,u,step_vec);
u2 = method(h_old/2,i-1,u,step_vec); % extrapolation to estimate the error
if abs(u2)>10^200 | strcmp(num2str(u2),'NaN')
errordlg('Your method produces a solution which can not be displayed. Change the chosen parameters!','error')
stop_calculation=1;
break
end
est=abs((2^p*u2-u1)/(2^p-1)-u1);
if est<TOL & est>0.1*TOL % step size will be accepted
RESULT.steps(i)=h_old;
u(i)=u1;
i=i+1;
else
h_new=h_old*(0.5*TOL/est)^(1/(p+1));
if h_new<h_old % decrease of the step size
h_old=h_new;
elseif h_new>=h_old & est<0.1*TOL % increase of the step size for the next step
u(i)=u1;
RESULT.steps(i)=h_old;
h_old=h_new;
i=i+1;
end
end
if step_vec(i-1)>SETTINGS.endvalue % cancel because the interval end is crossed
go = false;
end
end
else
start=i;
for i=start:length(step_vec)
u(i)=method(SETTINGS.step,i-1,u,step_vec);
if abs(u(i))>10^200 | strcmp(num2str(u(i)),'NaN')
errordlg('Your method produces a solution which can not be displayed. Change the chosen parameters!','error')
stop_calculation=1;
break
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%% Die Stabilittsgebiete %%%%%%%%%%%%%%%%%%%%%%%
if SETTINGS.calculate_region_of_stability==1
% der Bereich der auf Stabilitt berprft wird
X = -5:0.05:2;
Y = -4:0.05:4;
i = sqrt(-1);
count_idx = 0;
if SETTINGS.isRKV
% Runge-Kutta-verfahren
for x_idx = 1:length(X)
for y_idx = 1:length(Y)
if abs(1+c*(eye(s)-(X(x_idx)+i*Y(y_idx))*B)^(-1)*ones(1,s)'*(X(x_idx)+i*Y(y_idx)))<=1 % die Stabilitts Funktion eines RKV in allgemeiner Form wird
count_idx = count_idx+1; % fr die verschiedenen Gitterpunkte ausgewertet.
plot_vec_x(count_idx) = X(x_idx);
plot_vec_y(count_idx) = Y(y_idx);
end
end
end
else
% linear multi step methods
for x_idx = 1:length(X)
for y_idx = 1:length(Y)
Wurzeln = roots(rho-(X(x_idx)+i*Y(y_idx))*sigma);
if abs(Wurzeln) <= 1 & (length(find(Wurzeln==1))==1 | length(find(Wurzeln==1))==0)% Wurzelbedingung prfen
count_idx = count_idx+1;
plot_vec_x(count_idx) = X(x_idx);
plot_vec_y(count_idx) = Y(y_idx);
end
end
end
end
if exist('plot_vec_y','var')
RESULT.x = plot_vec_x;
RESULT.y = plot_vec_y;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if stop_calculation == 1
RESULT.numerical_solution = 0;
RESULT.gridpoints = 0;
RESULT.exact_solution = 0;
RESULT.steps = 0;
else
RESULT.numerical_solution = u;
if isfield(SETTINGS,'fu_exact')
% exact solution is known:
fu_exact = SETTINGS.fu_exact;
for j=1:length(step_vec)
u_exact(j) = fu_exact(step_vec(j));
end
RESULT.exact_solution = u_exact;
end
RESULT.gridpoints = step_vec;
if SETTINGS.isRKV
SETTINGS.B = B;
SETTINGS.c = c;
SETTINGS.a = a;
else
SETTINGS.rho = rho;
SETTINGS.sigma = sigma;
end
end
end
