ode45解微分方程（组）

87 阅读 0 评论 58 点赞

我是靠谱客的博主大胆芹菜，这篇文章主要介绍ode45解微分方程（组），现在分享给大家，希望可以做个参考。

适用非刚性方程

[T,Y] = ode45(odefun,tspan,y0)
解算 y′=f(t,y) 形式的方程组
[t,y] = ode45(odefun,tspan,y0,options) 还使用由 options（使用 odeset 函数创建的参数）定义的积分设置。例如，使用 AbsTol 和 RelTol 选项指定绝对误差容限和相对误差容限

odefun 是函数句柄，可以是函数文件名，或匿名函数句柄
tspan是区间[t0 tf] 或者一系列时间散点[t0,t1,…,tf]
y0是初始值向量
T返回列向量的时间点
Y返回对应T的求解列向量
生成的输出即为时间点 t 的列向量和解数组 y。y 中的每一行都与 t 的相应行中返回的时间相对应。y 的第一列与 $y_1$ 相对应，第二列与 $y_2$ 相对应。

ode45 仅适用于使用两个输入参数（t 和 y）的函数。但是，可以通过在外部定义参数并在指定函数句柄时传递这些额外参数。

复制代码

1
2
3
4
5
6
7
A = 1;
B = 2;
tspan = [0 5];
y0 = [0 0.01];
[t,y] = ode45(@(t,y)[y(2);(A/B)*t.*y(1)], tspan, y0);
plot(t,y(:,1),'-o',t,y(:,2),'-.')

- [example 1]

复制代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
function Testode45
tspan=[3.9 4.0]; %求解区间
y0=[8 2]; %初值
[t,y]=ode45(@odefun,tspan,y0);
plot(t,y(:,1),'-o',t,y(:,2),'-*')
legend('y1','y2')
title('y"=-t*y + e^t*y"+3sin2t')
xlabel('t')
ylabel('y')
end
function dx=odefun(t,x)
dx=zeros(2,1); % 列向量
dx(1)=x(2);
dx(2)=-t*x(1)+exp(t)*x(2)+3*sin(2*t); %常微分方程公式
end

或

复制代码

1
2
3
4
5
6
7
8
9
tspan=[3.9 4.0]; %求解区间
y0=[8 2]; %初值
[t,y]=ode45(@(t,y)[y(2);-t*y(1)+exp(t)*y(2)+3*sin(2*t)],tspan,y0);
plot(t,y(:,1),'-o',t,y(:,2),'-*')
legend('y1','y2')
title('y"=-t*y + e^t*y"+3sin2t')
xlabel('t')
ylabel('y')

或

复制代码

1
2
3
4
5
6
7
8
9
10
tspan=[3.9 4.0]; %求解区间
y0=[8 2]; %初值
fcn=@(t,y)[y(2);-t*y(1)+exp(t)*y(2)+3*sin(2*t)];
[t,y]=ode45(fcn,tspan,y0);
plot(t,y(:,1),'-o',t,y(:,2),'-*')
legend('y1','y2')
title('y"=-t*y + e^t*y"+3sin2t')
xlabel('t')
ylabel('y')

- [ example 2]

复制代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
% 表示该方程组
m = 1;
dy = @(t, y)[y(2); (5 - y(2)*exp(y(1)) + y(1)^2)/m];
y0=[1;-10];
tspan = [0, 5];
% m = 1,质量为1
[t_m_1, y_m_1] = ode45(dy, tspan, y0);
% m = 2质量为2
m = 2;
dy = @(t, y)[y(2); (5 - y(2)*exp(y(1)) + y(1)^2)/m];
[t_m_2, y_m_2] = ode45(dy, tspan, y0);
% plot
subplot(1, 2, 1);
plot(t_m_1, y_m_1(:, 2));
hold on
plot(t_m_2, y_m_2(:, 2));
title('dy/dt')
legend('m=1','m=2')
subplot(1, 2, 2);
plot(t_m_1, y_m_1(:, 1));
hold on
plot(t_m_2, y_m_2(:, 1));
title('y')
legend('m=1','m=2')

复制代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
m = 1;
y0=[1;-10];
tspan = [0, 5];
[t_m_1, y_m_1] = ode45(@(t, y)[y(2); (5 - y(2)*exp(y(1)) + y(1)^2)/m], tspan, y0);
% m = 2
m = 2;
[t_m_2, y_m_2] = ode45(@(t, y)[y(2); (5 - y(2)*exp(y(1)) + y(1)^2)/m], tspan, y0);
% plot
subplot(1, 2, 1);
plot(t_m_1, y_m_1(:, 2));
hold on
plot(t_m_2, y_m_2(:, 2));
title('dy/dt')
legend('m=1','m=2')
subplot(1, 2, 2);
plot(t_m_1, y_m_1(:, 1));
hold on
plot(t_m_2, y_m_2(:, 1));
title('y')
legend('m=1','m=2')