我是靠谱客的博主 大胆芹菜,最近开发中收集的这篇文章主要介绍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')

最后

以上就是大胆芹菜为你收集整理的ode45解微分方程(组)的全部内容,希望文章能够帮你解决ode45解微分方程(组)所遇到的程序开发问题。

如果觉得靠谱客网站的内容还不错,欢迎将靠谱客网站推荐给程序员好友。

本图文内容来源于网友提供,作为学习参考使用,或来自网络收集整理,版权属于原作者所有。
点赞(41)

评论列表共有 0 条评论

立即
投稿
返回
顶部