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_1 y1 相对应,第二列与 y 2 y_2 y2 相对应。

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

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]

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

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')

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]

% 表示该方程组
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')
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')
  • 6
    点赞
  • 85
    收藏
    觉得还不错? 一键收藏
  • 2
    评论
评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值