MATLAB ODE start / stop conditions for ode15i

I was looking for ways to stop MATLAB ode when it meets a certain condition. I found the answer in this section of the MatLab ODE Start / Stop Conditions , where the use of "events" was discussed. However, this applies to ode45, and when I tried to use the โ€œeventsโ€ with ode15i, it just didn't work, and MATLAB showed an error.

I am trying to learn this with a simple example and solved a simple system of differential equations, as shown below.

dx / dt = 5x + 3y; dy / dt = x + 7y; I decided to use them ode45 and tried to do the same with ode15i, but this will not work. Below are my codes.

With ode45

    function start_stop_test_ode45
    y0 = [5;1];
    tv = linspace(0,2,100);
    options = odeset('Events',@events);    
    f = @(t,y) [5*y(1) + 3*y(2);y(1) + 7*y(2)];
    [t,Y] = ode45(f,tv,y0,options);
    xNI = Y(:,1);
    yNI = Y(:,2);
    xCF   =  3*exp(4*t) + 2*exp(8*t);
    yCF   = -1*exp(4*t) + 2*exp(8*t);
    % Here we plot all the graphs
    figure(1)
    plot(t,xNI,'--k',t,xCF,'r','Linewidth',1.75)
    xlabel('t (s)')
    ylabel('x')
    legend('Numerical Solution','Closed Form Solution')
    figure(2)
    plot(t,yNI,'--k',t,yCF,'r','Linewidth',1.75)
    xlabel('t (s)')
    ylabel('y')
    legend('Numerical SOlution','Closed Form Solution')
    % Here we solve plot the variation of x with y
    figure(3)
    plot(xNI,yNI,'k','Linewidth',2);
    end

    function [value,isterminal,direction] = events(t,y)
    value = [y(1) - 7782;y(2) - 8863];  % Detect y = 7356
    isterminal = [1;1];
    direction = [0;0];
    end

With ode15i

    function start_stop_test_ode15i
    clc;clear all
    t0 = 0;
    y0 = [5;1];
    Fxdy0 = [1;1];
    Fxdyp0 = [0;0];
    yp0 = [28;12];
    tRange = [0 2];
    options = odeset('Events',@events);
    [y0,yp0] = decic(@ode15ifun,t0,y0,Fxdy0,yp0,Fxdyp0);
    sol = ode15i(@ode15ifun,tRange,y0,yp0,options);
    tv = linspace(0,2,100);
    sv = deval(sol,tv);
    sv = sv';
    t = tv;
    xNI = sv(:,1);
    yNI = sv(:,2);
    xCF   =  3*exp(4*t) + 2*exp(8*t);
    yCF   = -1*exp(4*t) + 2*exp(8*t);
    % Here we plot all the graphs
    figure(4)
    plot(t,xNI,'--k',t,xCF,'r','Linewidth',1.75)
    xlabel('t (s)')
    ylabel('x')
    legend('Numerical Solution','Closed Form Solution')
    figure(5)
    plot(t,yNI,'--k',t,yCF,'r','Linewidth',1.75)
    xlabel('t (s)')
    ylabel('y')
    legend('Numerical SOlution','Closed Form Solution')
    % Here we solve plot the variation of x with y
    figure(6)
    plot(xNI,yNI,'k','Linewidth',2);
    end

    function [value,isterminal,direction] = events(t,y)
    value = [y(1) - 7782;y(2) - 8863];  % Detect y = 7356
    isterminal = [1;1];
    direction = [0;0];
    end

Where is ode15ifun

    function res = ode15ifun(t,y,yp)
    %UNTITLED3 Summary of this function goes here
    %   Detailed explanation goes here
        res1 = yp(1) - 5*y(1)- 3*y(2);
        res2 = yp(2) - y(1) - 7*y(2);
        res  = [res1;res2];
    end

ode45 , ode15i . - , ode15i?

( TroyHaskin) . enter image description here

+4

Source: https://habr.com/ru/post/1607827/


All Articles