Non-stiff ODE在MATLAB中非刚性ODE求解

 1 ode45  Solve non-stiff differential equations, medium order method.
 2  [TOUT,YOUT] = ode45(ODEFUN,TSPAN,Y0) integrates the system of
 3  differential equations y' = f(t,y) from time TSPAN(1) to TSPAN(end)
 4  with initial conditions Y0. Each row in the solution array YOUT
 5  corresponds to a time in the column vector TOUT. 
 6    * ODEFUN is a function handle. For a scalar T and a vector Y,
 7    ODEFUN(T,Y) must return a column vector corresponding to f(t,y).
 8    * TSPAN is a two-element vector [T0 TFINAL] or a vector with
 9    several time points [T0 T1 ... TFINAL]. If you specify more than
10    two time points, ode45 returns interpolated solutions at the
11    requested times.
12    * YO is a column vector of initial conditions, one for each equation.
13 
14  [TOUT,YOUT] = ode45(ODEFUN,TSPAN,Y0,OPTIONS) specifies integration
15  option values in the fields of a structure, OPTIONS. Create the
16  options structure with odeset.
17 
18  [TOUT,YOUT,TE,YE,IE] = ode45(ODEFUN,TSPAN,Y0,OPTIONS) produces
19  additional outputs for events. An event occurs when a specified function
20  of T and Y is equal to zero. See ODE Event Location for details.
21 
22  SOL = ode45(...) returns a solution structure instead of numeric
23  vectors. Use SOL as an input to DEVAL to evaluate the solution at
24  specific points. Use it as an input to ODEXTEND to extend the
25  integration interval.
26 
27  ode45 can solve problems M(t,y)*y' = f(t,y) with mass matrix M that is
28  nonsingular. Use ODESET to set the 'Mass' property to a function handle
29  or the value of the mass matrix. ODE15S and ODE23T can solve problems
30  with singular mass matrices.
31 
32  ODE23, ode45, ODE78, and ODE89 are all single-step solvers that use
33  explicit Runge-Kutta formulas of different orders to estimate the error
34  in each step.
35    * ode45 is for general use.
36    * ODE23 is useful for moderately stiff problems.
37    * ODE78 and ODE89 may be more efficient than ode45 on non-stiff problems
38    that are smooth except possibly for a few isolated discontinuities.
39    * ODE89 may be more efficient than ODE78 on very smooth problems, when 
40    integrating over long time intervals, or when tolerances are tight.
41 
42  Example
43      [t,y]=ode45(@vdp1,[0 20],[2 0]);   
44      plot(t,y(:,1));
45    solves the system y' = vdp1(t,y), using the default relative error
46    tolerance 1e-3 and the default absolute tolerance of 1e-6 for each
47    component, and plots the first component of the solution. 
48  
49  Class support for inputs TSPAN, Y0, and the result of ODEFUN(T,Y):
50    float: double, single
51 
52  See also ode23, ode78, ode89, ode113, ode15s, ode23s, ode23t, ode23tb,
53       ode15i, odeset, odeplot, odephas2, odephas3, odeprint, deval,
54       odeexamples, function_handle.
55
56  Documentation for ode45
57  Other uses of ode45

1[t, y] = ode45(@odefun, tspan, y0);

1sol = ode45(@odefun, tspan, y0);

1ode45(@odefun, tspan, y0);

1options = odeset('RelTol',1e-6,'AbsTol',1e-9);
2[t, y] = ode45(@odefun, tspan, y0, options);

1[t, y] = ode45(@(t, y) odefun(t, y, param), tspan, y0);

1[t, y] = ode45(@odefun, tspan, y0, [], param);

1function dydt = odefun(t, y, param)
2  % 这里的param可以是一个自定义参数
3  dydt = zeros(2, 1); % 假设有两个变量
4  dydt(1) = y(2); % y1' = y2
5  dydt(2) = -param * y(1); % y2' = -param * y1
6end

1function dydt = rayleigh_ode(t, y, mu)
2% Rayleigh ODE function
3% dydt = rayleigh(t, y, mu)
4% t: time variable
5% y: state variable (y(1) = y, y(2) = dy/dt)
6dydt = zeros(2, 1);
7dydt(1) = y(2);
8dydt(2) = mu * (1 - (1/3) * y(2)^2) * y(2) - y(1);
9end

1function dydt = odefun(t, y, mu)
2dydt = [0, 1;...
3  -1, mu * (1 - (1/3) * y(2)^2)] * y(:);
4end

1odefun = @(t, y, mu)[0, 1;...
2  -1, mu * (1 - (1/3) * y(2)^2)] * y(:);
3
4ode45(odefun, [0 100], [0.1; 0], [], 0.1);

1arguments
2    mu (1,1) double {mustBeNumeric, mustBeFinite} = 0.4; % Default value for mu
3    show_fig (1,1) logical = false; % Default value for show_fig
4    y_1_0 (1,1) double {mustBeNumeric, mustBeFinite} = 0.001;
5end

1tl = tiledlayout(1,2);
2title(tl, sprintf('Rayleigh ODE Solution mu = %.2f', mu));
3
4nexttile
5plot(t, y(:, 1))
6xlabel('Time(s)')
7ylabel('y') 

 1function rayleigh(mu, show_fig, y_1_0)
 2arguments
 3    mu (1,1) double {mustBeNumeric, mustBeFinite} = 0.4; % Default value for mu
 4    show_fig (1,1) logical = false; % Default value for show_fig
 5    y_1_0 (1,1) double {mustBeNumeric, mustBeFinite} = 0.001;
 6end
 7
 8y0 = [y_1_0; 0]; % Initial conditions: y(0) = 0, dy/dt(0) = 1
 9tspan = [0 100]; % Time span for the solution
10
11
12fig = figure('Units','inches', 'Position', [1, 1, 12, 6], 'Visible', false);
13
14% [t, y] = ode45(@(t, y) rayleigh_ode(t, y, mu), tspan, y0);
15[t, y] = ode45(@rayleigh_ode, tspan, y0, [], mu);
16
17tl = tiledlayout(1,2);
18title(tl, sprintf('Rayleigh ODE Solution mu = %.2f', mu));
19
20nexttile
21plot(t, y(:, 1))
22xlabel('Time(s)')
23ylabel('y')
24
25nexttile
26plot(y(:,1), y(:,2), 'b-', 'LineWidth', 2);
27xlabel('y')
28ylabel('y''')
29
30exportgraphics(tl, sprintf('rayleigh_solution_%.2f.png', mu), 'Resolution', 500);
31
32if show_fig
33    set(fig, 'Visible', 'on');
34else
35    close(fig); % Close the figure if not showing
36end
37end
38
39function dydt = rayleigh_ode(t, y, mu)
40% Rayleigh ODE function
41% dydt = rayleigh(t, y, mu)
42% t: time variable
43% y: state variable (y(1) = y, y(2) = dy/dt)
44dydt = zeros(2, 1);
45dydt(1) = y(2);
46dydt(2) = mu * (1 - (1/3) * y(2)^2) * y(2) - y(1);
47end