Why ODE15s does not work like an Euler's method with fixed point?

5 次查看（过去 30 天）

Cesar García Echeverry 2020-4-26

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/520902-why-ode15s-does-not-work-like-an-euler-s-method-with-fixed-point

评论： Cesar García Echeverry 2020-4-26

When i run this code for Euler´s method, it runs Ok, but when i try to use ODE15s it shows a wrong answer. It is ODE15s well defined?. I'll appreciate your help.

function []=cr_NRTL(X0,x1,x2,l,L1)
global X1 X2 Xea estado
x1inicial=x1;
x2inicial=x2;
sol_daes=2; %1 Euler, 2 ODE
tspan=[0 l];
opts=odeset('NormControl','on','AbsTol',1e-10,'RelTol',1e-5);
for k=1:2
    x1=x1inicial;
    x2=x2inicial;
    for i=1:L1
        X1 = x1; X2 = x2;
        
        %Para solucionar la envolvente
        options = optimset('Display','off');
        lb=[0,0,0,0,0,0,0,0,300];
        ub=[1,1,1,1,1,1,1,1,400];
        Xea=lsqnonlin(@funcion,X0,lb,ub,options);
        
        X0=Xea;
        
        for j=1:9
            xea(i,j)=Xea(j);
        end
        
        %Para evaluar las residuales del ELV fuera de la envolvente
        if k==1 && (Xea(7)>0.9999 || Xea(8)>0.9999)
            if Xea(7)>1 || Xea(8)>1
                Xea
            end
            cr_NRTL2(x1,x2,l);
            break
        else
            
            xglobala=[x1; x2];
            
            %Para solucionar las residuales del ELLV
            switch sol_daes
                
                case 1 %Euler
                    switch k
                        case 1
                            tao=0.02;
                            dxdt=cambioELLV(xglobala);
                            xsa=[xglobala(1)+dxdt(1)*tao;
                                xglobala(2)+dxdt(2)*tao];
                        case 2
                            tao=0.1;
                            dxdt=cambioELLV(xglobala);
                            xsa=[xglobala(1)-dxdt(1)*tao;
                                xglobala(2)-dxdt(2)*tao];
                    end
                    
                case 2 %ODE
                    switch k
                        case 1
                            [~,dxdt]=ode15s(@cambioELLV,tspan,xglobala,opts);
                        case 2
                            [~,dxdt]=ode15s(@cambioELLV,-tspan,xglobala,opts);
                    end
                    xsa=dxdt;
                    
            end
            
            x1=xsa(1);
            x2=xsa(2);
            
            %                     if k==2
            for j=1:2
                xsalida(i,j)=xglobala(j);
                %                             end
            end
            
        end
        
        %             end
        %             if k ==2
        xe = xea;
        %             end
    end
    toc
    % Grafica
    hold on
    if estado==1
        plot(xe(:,1),xe(:,2),'r')
    end
    plot(xsalida(:,1),xsalida(:,2),'b');
end
end