Description of ivp

0001 function [RESULT, SETTINGS]=ivp_solve(SETTINGS)
0002 % function thal solves the initial value problem
0003 % by using the specified parameters out of SETTINGS
0004 %
0005 % input:   o SETTINGS....struct containing all information about the
0006 %                        current program run
0007 %
0008 % example  o SETTINGS =
0009 %
0010 %                              step: 0.1000
0011 %                            method: {'explicit Euler-method'}
0012 %     calculate_region_of_stability: 0
0013 %                          equation: @(x,y)(cos(y)).^2./sin(y).*x.*exp(1/2.*x)
0014 %                          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)))
0015 %                        beginvalue: -5
0016 %                          endvalue: 4
0017 %                        startvalue: 1.3082
0018 %                             isRKV: 1
0019 %                      step_control: 0
0020 %
0021 % output:  o RESULT....struct containing all needed information about the
0022 %                      numerical soltuion
0023 %
0024 % example: o RESULT =
0025 %
0026 %     numerical_solution: [1x91 double]
0027 %         exact_solution: [1x91 double]
0028 %             gridpoints: [1x91 double]
0029 %
0030 % notes:   o
0031 %
0032 % authors: o Christian Jäkel (University of Technology Dresden)
0033 %          o Daniel Klawitter (University of Technology Dresden)
0034 
0035 
0036 stop_calculation = 0;
0037 f = SETTINGS.equation;
0038 
0039 % calculates the k's for implicit RKV
0040     function [ka] = k_imp(j,age)
0041         function[Funktion] = F(x,y1,y2)
0042             Funktion       = y2-f(x+a*age,(y1'+age*B*y2')');
0043         end
0044         ka = fsolve(@(y)F(step_vec(j),u(j),y),ones(1,s),optimset('Display','off'));
0045         
0046     end
0047 
0048 
0049 % calculates the k's for explicit RKV
0050     function [ka] = k(j,age)
0051         for i2 = 1:s
0052             ka(1)  = f(step_vec(j)+a(i2)*age,u(j));
0053             ka(i2) = f(step_vec(j)+a(i2)*age,u(j)+age*sum(B(i2,1:i2-1).*ka(1:i2-1)));
0054         end
0055     end
0056 
0057 % this function calculates initial values for the linear multistep methods
0058     function[startwerte] = berechne_Startwerte(stufe)                                           
0059         a = [0, 1/5, 3/10, 4/5, 8/9, 1, 1];
0060         B = [0        ,           0,          0,        0,           0,     0,   0;
0061             1/5       ,           0,          0,        0,           0,     0,   0;
0062             3/40      ,        9/40,          0,        0,           0,     0,   0;
0063             44/45     ,      -56/15,       32/9,        0,           0,     0,   0;
0064             19372/6561, -25360/2187, 64448/6561, -212/729,           0,     0,   0;
0065             9017/3168 ,     -355/33, 46732/5247,   49/176, -5103/18656,     0,   0;
0066             35/384    ,           0,   500/1113,  125/192,  -2187/6784, 11/84,   0];
0067         c = [5179/57600,0,7571/16695,393/640,-92097/339200,187/2100,1/40];
0068         s = 7;
0069         p=5;
0070         tele(1)=SETTINGS.startvalue;
0071         
0072         function [ka] = kaa(j,age)
0073             for i2 = 1:s
0074                 ka(1)  = f(step_vec(j)+a(i2)*age,tele(j));
0075                 ka(i2) = f(step_vec(j)+a(i2)*age,tele(j)+age*sum(B(i2,1:i2-1).*ka(1:i2-1)));
0076             end
0077         end
0078         if SETTINGS.step_control==1
0079             i=2;
0080             method=@(sta,z,vector)(vector(z)+sta*sum(c.*kaa(z,sta)));
0081             while go==true
0082                 u1 = method(h_old,i-1,tele);
0083                 u2 = method(h_old/2,i-1,tele);  % extrapolation to estimate the error
0084                 est=abs((2^p*u2-u1)/(2^p-1)-u1);
0085                 if est<TOL & est>0.1*TOL   % step size will be accepted
0086                     RESULT.steps(i)=h_old;
0087                     step_vec(i)=step_vec(i-1)+h_old;
0088                     tele(i)=u1;
0089                     i=i+1;
0090                 else
0091                     h_new=h_old*(0.5*TOL/est)^(1/(p+1));
0092                     if h_new<h_old  % decrease of the step size
0093                         h_old=h_new;
0094                         
0095                     elseif h_new>=h_old & est<0.1*TOL    % increase of the step size for the next step
0096                         tele(i)=u1;
0097                         step_vec(i)=step_vec(i-1)+h_old;
0098                         RESULT.steps(i)=h_old;
0099                         h_old=h_new;
0100                         i=i+1;
0101                     end
0102                 end
0103                 if i==length(rho)
0104                     go = false;
0105                     
0106                 end
0107                 
0108             end
0109         else
0110             for i = 2:stufe
0111                 tele(i) = tele(i-1)+SETTINGS.step*sum(c.*kaa(i-1,SETTINGS.step));
0112             end
0113         end
0114         
0115         startwerte = tele;
0116     end
0117 
0118 
0119 
0120 %%%%%%%%%%%%%%%%%%%%%%%%%%%%     THE METHODS    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0121 switch char(SETTINGS.method)
0122     case {'explicit Euler-method'}
0123         % Euler explizit  --> explizit
0124         B = [0];
0125         a = [0];
0126         c = [1];
0127         s = 1;
0128         p=1;
0129     case {'Heun-Verfahren'}
0130         % Heun-Verfahren --> explizit
0131         B = [0, 0;
0132             1, 0];
0133         a = [0  , 1  ];
0134         c = [1/2, 1/2];
0135         s = 2;
0136         p=2;
0137     case {'Klassisches RKV'}
0138         % klassisches RKV --> explizit
0139         B = [0 , 0 , 0, 0;
0140             1/2, 0 , 0, 0;
0141             0  ,1/2, 0, 0;
0142             0  ,  0, 1, 0];
0143         a = [0  , 1/2, 1/2 ,  1];
0144         c = [1/6, 1/3, 1/3, 1/6];
0145         s = 4;
0146         p=4;
0147     case {'Dormand-Prince RKV'}
0148         % embedded RKV from Dormand-Prince sheme
0149         a = [0, 1/5, 3/10 , 4/5, 8/9, 1 , 1 ];
0150         B = [0        , 0          , 0          , 0       , 0            , 0       , 0;
0151             1/5       , 0          , 0          , 0       , 0            , 0       , 0;
0152             3/40      , 9/40       , 0          , 0       , 0            , 0       , 0;
0153             44/45     , -56/15     , 32/9       , 0       , 0            , 0       , 0;
0154             19372/6561, -25360/2187, 64448/6561 , -212/729, 0            , 0       , 0;
0155             9017/3168 , -355/33    , 46732/5247 , 49/176  , -5103/18656  , 0       , 0;
0156             35/384    , 0          , 500/1113   , 125/192 , -2187/6784   , 11/84   , 0];
0157         % with this c order4
0158         c = [5179/57600 , 0, 7571/16695 , 393/640, -92097/339200, 187/2100, 1/40];
0159         % and with this one order 5
0160         c = [35/384    , 0          , 500/1113   , 125/192 , -2187/6784   , 11/84   , 0];
0161         s = 7;
0162         p=5;
0163     case {'implicit Euler-method'}
0164         % implizites verfahren
0165         B = [1];
0166         a = [1];
0167         c = [1];
0168         s = 1;
0169         p=1;
0170     case {'Radau II RKV'}
0171         % Radau IIA method order 3 entnommen von
0172         % http://www-user.tu-chemnitz.de/~benner/Lehre/NumerikODE/Folie_IRK.pdf
0173         B = [5/12, -1/12;
0174             3/4 ,  1/4 ];
0175         a = [1/3,1];
0176         c = [3/4,1/4];
0177         s = 2;
0178         p=3;
0179     case {'Gauß RKV'}
0180         % gaußverfahren der ordnung 6 entnommen von
0181         % http://www-user.tu-chemnitz.de/~benner/Lehre/NumerikODE/Folie_IRK.pdf
0182         B = [5/36              , (10-3*sqrt(15))/45, (25-6*sqrt(15))/180;
0183             (10+3*sqrt(15))/72 ,  2/9              , (10-3*sqrt(15))/72;
0184             (25+6*sqrt(15))/180, (10+3*sqrt(15))/45,  5/36];
0185         a = [(5-sqrt(15))/10,1/2,(5+sqrt(15))/10];
0186         c = [5/18,4/9,5/18];
0187         s = 3;
0188         p=6;
0189 end
0190 %%%%%%%%%%%%%%%%%     calculation of the numerical solution      %%%%%%%%%%%%%%%%%%%%%%
0191 
0192 step_vec(1)=SETTINGS.beginvalue;  % vector with grid points
0193 h_old=SETTINGS.step;
0194 if SETTINGS.step_control==1
0195     RESULT.steps(1)=SETTINGS.step;  % vektor with step size for each grid point
0196     go=true;
0197     TOL=SETTINGS.tolerance;
0198 else
0199     step_vec = SETTINGS.beginvalue:SETTINGS.step:SETTINGS.endvalue;
0200 end
0201 
0202 switch char(SETTINGS.method)
0203     case {'explicit Euler-method','Heun-Verfahren','Klassisches RKV','Dormand-Prince RKV'}
0204         % explizite RKV
0205         u(1) = SETTINGS.startvalue;
0206         method=@(sta,z,vector,vector_step)(vector(z)+sta*sum(c.*k(z,sta)));
0207         i=2;
0208     case {'implicit Euler-method','Radau II RKV','Gauß RKV','radauIIA'}
0209         u(1) = SETTINGS.startvalue;
0210         method=@(sta,z,vector,vector_step)(vector(z)+sta*sum(c.*k_imp(z,sta)));
0211         i=2;
0212     case {'Adams-Bashforth(2)'}
0213         % ein zweistufiges Adamsverfahren (explizit)
0214         rho        = [1, -1 ,  0  ];
0215         sigma      = [0, 3/2, -1/2];
0216         u(1)       = SETTINGS.startvalue;
0217         startwerte = berechne_Startwerte(2);
0218         u(2)       = startwerte(2);
0219         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))));
0220         i=3;
0221         go=true;
0222     case {'Adams-Bashforth(4)'}
0223         % ein vierstufiges Verfahren (explizit)
0224         rho        = [1, -1   , 0     , 0 ,0  ];
0225         sigma      = [0, 55/24, -59/24, 37/24,-9/24];
0226         u(1)       = SETTINGS.startvalue;
0227         startwerte = berechne_Startwerte(4);
0228         u(2:4)     = startwerte(2:4);
0229         i=5;
0230         go=true;
0231         method = @(sta,z,vector,vector_step)(vector(z)+sta/24*(55*f(vector_step(z),vector(z))...
0232             -59*f(vector_step(z-1),vector(z-1))...
0233             +37*f(vector_step(z-2),vector(z-2))...
0234             -9*f(vector_step(z-3),vector(z-3))));
0235         
0236     case {'Adams-Moulton(2)'}
0237         % Adams-Moulton-verfahren, implizit, 2-stufig
0238         rho        = [1,-1,0];
0239         sigma      = [5/12,8/12,-1/12];
0240         u(1)       = SETTINGS.startvalue;
0241         startwerte = berechne_Startwerte(2);
0242         u(2)       = startwerte(2);
0243         i=3;
0244         go=true;
0245         method = @(sta,z,vector,vector_step)(fzero(@(y) y-vector(z)-sta/12*(5*f(vector_step(z+1),y)...
0246             +8*f(vector_step(z),vector(z))...
0247             -f(vector_step(z-1),vector(z-1))),vector(z)));
0248         
0249     case {'Milne-Simpson(4)'}
0250         % Milne-Simpson-verfahren, implizit, 4-stufig
0251         rho        = [1    , 0     , -1   , 0   , 0    ];
0252         sigma      = [29/90, 124/90, 24/90, 4/90, -1/90];
0253         u(1)       = SETTINGS.startvalue;
0254         startwerte = berechne_Startwerte(4);
0255         u(2:4)     = startwerte(2:4);
0256         i=5;
0257         go=true;
0258         method = @(sta,z,vector,vector_step) (fzero(@(y)y-vector(z-1)-sta/90*(29*f(vector_step(z+1),y)+...
0259             124*f(vector_step(z),vector(z))...
0260             +24*f(vector_step(z-1),vector(z-1))...
0261             +4*f(vector_step(z-2),vector(z-2))...
0262             -f(vector_step(z-3),vector(z-3))),vector(z)));
0263         
0264     case {'BDF 5'}
0265         % BDF-Verfahren, implizit, 5-stufig
0266         rho        = [137/60, -300/60, 300/60, -200/60, 75/60, -12/60];
0267         sigma      = [1     , 0      , 0     , 0      , 0    ,      0];
0268         u(1)       = SETTINGS.startvalue;
0269         startwerte = berechne_Startwerte(5);
0270         u(2:5)     = startwerte(2:5);
0271         i=6;
0272         go=true;
0273         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))...
0274             -sta*f(vector_step(z+1),y),vector(z)));
0275         
0276     case {'BDF 2'}
0277         % BDF Verfahren 2-stufig, A-Stabil!!!!
0278         rho=[3/2,-4/2,1/2];
0279         sigma=[1,0,0];
0280         u(1) = SETTINGS.startvalue;
0281         startwerte = berechne_Startwerte(2);
0282         u(2)=startwerte(2);
0283         i=3;
0284         go=true;
0285         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)));
0286         
0287     case'eigenesRKV'
0288         
0289         B = SETTINGS.B;
0290         c = SETTINGS.c;
0291         a = SETTINGS.a;
0292         s = length(B);
0293         u(1) = SETTINGS.startvalue;
0294         i=2;
0295         p=length(B); % rough estimation for the order of convergence.
0296         method=@(sta,z,vector,vector_step)(vector(z)+sta*sum(c.*k_imp(z,sta)));
0297         
0298         
0299         %if abs(u(i))>10^200
0300         %      errordlg('Your method produces a solution which can not be displayed. Change the chosen parameters!','error')
0301         %      stop_calculation=1;
0302         %      break
0303         %
0304         %  end
0305         
0306     case'eigenesMSV'
0307         
0308         u(1)   = SETTINGS.startvalue;
0309         rho    = SETTINGS.rho;
0310         sigma  = SETTINGS.sigma;
0311         rho2   = fliplr(rho);
0312         sigma2 = fliplr(sigma);
0313         st     = length(rho)-1;
0314         startwerte = berechne_Startwerte(st);
0315         u(2:st)= startwerte(2:st);
0316         i=st+1;
0317         go=true;
0318         p=length(rho);
0319         
0320         
0321         %zwischenwert1 = sum(rho2(1:st).*u(i-st:i-1));
0322         %zwischenwert2 = sum(sigma2(1:st).*f(step_vec(i-st:i-1),u(i-st:i-1)));
0323         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)));
0324         
0325         %if abs(u(i))>10^200 | strcmp(num2str(u(i)),'NaN')
0326         %    errordlg('Your method produces a solution which can not be displayed. Change the chosen parameters!','error')
0327         %   stop_calculation=1;
0328         %   break
0329         %
0330         %end
0331         
0332         
0333 end
0334 
0335 if SETTINGS.step_control==1  % step control
0336     
0337     while go==true
0338         step_vec(i)=step_vec(i-1)+h_old; % this is for the implicit methods because they need this point for the error estimation
0339         u1 = method(h_old,i-1,u,step_vec);
0340         u2 = method(h_old/2,i-1,u,step_vec);  % extrapolation to estimate the error
0341         
0342         if abs(u2)>10^200 | strcmp(num2str(u2),'NaN')
0343             errordlg('Your method produces a solution which can not be displayed. Change the chosen parameters!','error')
0344             stop_calculation=1;
0345             break
0346             
0347         end
0348         
0349         est=abs((2^p*u2-u1)/(2^p-1)-u1);
0350         if est<TOL & est>0.1*TOL   % step size will be accepted
0351             RESULT.steps(i)=h_old;
0352             u(i)=u1;
0353             i=i+1;
0354         else
0355             h_new=h_old*(0.5*TOL/est)^(1/(p+1));
0356             if h_new<h_old  % decrease of the step size
0357                 h_old=h_new;
0358                 
0359             elseif h_new>=h_old & est<0.1*TOL    % increase of the step size for the next step
0360                 u(i)=u1;
0361                 RESULT.steps(i)=h_old;
0362                 h_old=h_new;
0363                 i=i+1;
0364             end
0365         end
0366         
0367         if step_vec(i-1)>SETTINGS.endvalue    % cancel because the interval end is crossed
0368             go = false;
0369             
0370         end
0371         
0372     end
0373     
0374 else
0375     
0376     start=i;
0377     for i=start:length(step_vec)
0378         u(i)=method(SETTINGS.step,i-1,u,step_vec);
0379         if abs(u(i))>10^200 | strcmp(num2str(u(i)),'NaN')
0380             errordlg('Your method produces a solution which can not be displayed. Change the chosen parameters!','error')
0381             stop_calculation=1;
0382             break
0383             
0384         end
0385     end
0386     
0387     
0388 end
0389 
0390 
0391 %%%%%%%%%%%%%%%%%%%%%%%%     Die Stabilitätsgebiete        %%%%%%%%%%%%%%%%%%%%%%%
0392 if SETTINGS.calculate_region_of_stability==1
0393     % der Bereich der auf Stabilität überprüft wird
0394     X = -5:0.05:2;
0395     Y = -4:0.05:4;
0396     i = sqrt(-1);
0397     
0398     count_idx = 0;
0399     if SETTINGS.isRKV
0400         % Runge-Kutta-verfahren
0401         for x_idx = 1:length(X)
0402             for y_idx = 1:length(Y)
0403                 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 Stabilitäts Funktion eines RKV in allgemeiner Form wird
0404                     count_idx             = count_idx+1;                                                          % für die verschiedenen Gitterpunkte ausgewertet.
0405                     plot_vec_x(count_idx) = X(x_idx);
0406                     plot_vec_y(count_idx) = Y(y_idx);
0407                 end
0408             end
0409         end
0410     else
0411         % linear multi step methods
0412         for x_idx = 1:length(X)
0413             for y_idx = 1:length(Y)
0414                 Wurzeln = roots(rho-(X(x_idx)+i*Y(y_idx))*sigma);
0415                 if abs(Wurzeln) <= 1 & (length(find(Wurzeln==1))==1 | length(find(Wurzeln==1))==0)% Wurzelbedingung prüfen
0416                     count_idx             = count_idx+1;
0417                     plot_vec_x(count_idx) = X(x_idx);
0418                     plot_vec_y(count_idx) = Y(y_idx);
0419                 end
0420             end
0421         end
0422     end
0423     
0424     if exist('plot_vec_y','var')
0425         RESULT.x = plot_vec_x;
0426         RESULT.y = plot_vec_y;  
0427     end
0428 end
0429 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0430 if stop_calculation == 1
0431     RESULT.numerical_solution = 0;
0432     RESULT.gridpoints = 0;
0433     RESULT.exact_solution = 0;
0434     RESULT.steps = 0;
0435 else
0436     RESULT.numerical_solution = u;
0437     if isfield(SETTINGS,'fu_exact')
0438         % exact solution is known:
0439         fu_exact = SETTINGS.fu_exact;
0440         for j=1:length(step_vec)
0441             u_exact(j) = fu_exact(step_vec(j));
0442         end
0443         RESULT.exact_solution = u_exact;
0444     end
0445     RESULT.gridpoints = step_vec;
0446     
0447     if SETTINGS.isRKV
0448         SETTINGS.B = B;
0449         SETTINGS.c = c;
0450         SETTINGS.a = a;
0451     else
0452         SETTINGS.rho = rho;
0453         SETTINGS.sigma = sigma;
0454     end
0455 end
0456 
0457 end
ivp_solve

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
