Matlab求解常微分方程组

头秃大嘴熊

已于 2022-12-21 10:44:41 修改

阅读量576

点赞数

文章标签： matlab 开发语言

于 2022-11-03 20:46:42 首次发布

本文链接：https://blog.youkuaiyun.com/cccczhh/article/details/124116768

版权

求解这个常微分方程组。

$\left\{ \begin{array}{l} \frac{d}{{dt}}(P{R^5}) = a{P^2}{R^5}\\ \\ \frac{{{d^2}R}}{{d{t^2}}} = {R^2P} \end{array} \right.$

初始条件为

$R(0) = \frac{1}{\varepsilon }$ $\frac { d R } { d t } ( 0 ) = - 1$ $\frac { d ^ { 2 } R } { d t ^ { 2 } } ( 0 ) = \varepsilon ^ { 3 }$ $P ( 0 ) = \varepsilon ^ { 5 } \$

其中ε取0.01，a是有上限的参数，求解方程的目的其实是找出a的临界值。

syms y(t) 
for i = [0:0.5:1.5,1.763]   %1.763是临界值
    eqn = diff(diff(y,2)*y^3) == i*diff(y,2)^2*y;


    %化三阶方程为一阶
    V = odeToVectorField(eqn);
    M = matlabFunction(V,'vars',{'t','Y'});
    e = 0.01;

    % 坐标范围和初始条件
    interval = [0 2/e];
    yInit = [1/e -1 e^3];
    ySol = ode45(M,interval,yInit);

    tValues =  linspace(0,2/e,1000);
    yValues = deval(ySol,tValues,1);
    pValues = deval(ySol,tValues,3)./yValues.^2;
    vValues = yValues.^3;
    plot(tValues,yValues)
    hold on
end

ylim([-10 100])
xlabel('$\hat{t}$','interpreter','latex');
ylabel('$\hat{R}$','interpreter','latex');
text(60,80,'\epsilon=0.01','FontName','Times New Roman','FontSize',15)
lgd = legend(num2str(0),num2str(0.5),num2str(1),num2str(1.5),num2str(1.763));
title(lgd,'\xi_n_o_\alpha')
set(gca,'LineWidth',2)
set(findall(gcf,'type','line'),'linewidth',2)
set(gca,'FontName','Times New Roman','FontSize',15)