时间戳

首页 > 解决方案 > 使用“subs”函数评估“dsolve”的输出在 Maltab 中提供额外的输出

matlab - 使用“subs”函数评估“dsolve”的输出在 Maltab 中提供额外的输出

问题描述

介绍

如果您想了解宏伟计划……请阅读简介。如果没有，请直接跳到我的问题。

我的微分方程和线性代数课程有一个项目，我必须使用计算机代数系统来求解具有常系数的一阶线性常微分方程和具有常系数的非线性 ODE。我必须证明这是在分析和数字上完成的。我的计划是有 2 个函数利用dsolve分析部分的函数和我通过 Matlab 的本教程ODE1了解的函数. 该函数被定义为使用欧拉方法逼近微分方程的解。这些函数都将使用相同的输入参数，因此每个输入都可以定义一次，并且这些函数都将了解要使用的输入（可能将函数嵌套在一个调用函数下）。目标是在用于查找数值近似值的相同区间内评估解析解，并在表格和图表中比较结果。我得到了数值解，以行向量的形式给了我一个“表格”，并绘制了该行向量。我开始遇到分析解决方案的问题...

我的问题

为了解决一阶线性 ODE，我生成了这个函数

function [s1] = L_Analytic(eqn,t0,h,numstep,y0)

% eqn is the differential equation to be solved
% t0 is the start of the interval that the resulting equation is to be evaluated at
% h is the stepsize
% numstep is the number of steps
% y0 is the initial condition

syms y(x)
cond = y(0) == y0;
A = dsolve(eqn,cond);
s1 = A;
S1 = s1;
   for t = t0 : h : h*(numstep-2)
       S1 = [subs(S1); vpa(subs(s1))]
   end 
end

这个函数生成的列表L_Analytic(diff(y)==y, 0, 0.1, 5, 1)是

1
1.0
1.105170...
1.221402...
1.349858...

在使用相同输入的 Matlab 中的不同函数中使用数值方法时，我得到列表：

1.0000
1.1000
1.2100
1.3310
1.4641

对于那些了解微分方程或精通微积分的人来说，y' = y 的解是 e^x，当使用 5 个步骤在 0:0.4 区间内进行评估时，列表应该是

1
1.105...
1.2214...
1.3498...
1.4918...

经过一些四舍五入。

1所以这里的问题是我的分析解决方案中有一个额外的条目。我相信它与循环中的subs(S1)部分有关S1 = [subs(S1); vpa(subs(s1))]，for但我不知道如何解决这个问题。

我有点理解为什么我需要使用该subs函数，因为我使用符号变量来使用dsolve在其答案中输出符号变量的函数。此外，为了for循环迭代和改变，符号变量必须替换为t每次的实际值。我确实尝试vpa(subs(s1))在循环之前移动for，但这只是在向量中返回了 5 次相同的值。我也试过不使用subs(S1)它给了我

exp(t)
1.0
1.1051709...
1.2214027...
1.3498588...

所以我很肯定这是代码的这一部分。

旁注：我了解分析方法输出列向量，如ODE1链接的视频中所示。为了让 Matlab 将其绘制为一条线，我将列向量转置为行向量，一旦解决方案部分固定，我将对解析解执行相同的操作。

标签： matlabmatrixfunctional-programmingdifferential-equationsdsolve

解决方案

通过改变for循环的内部结构，我让它工作了。我的最终函数代码是这样的：function [s1] = L_Analytic3(eqn,t0,h,numstep,y0)

%Differential Equation solver for specific inputs
%   eqn is the differential equation
%   t0 is start of evaluation interval
%   h is stepize
%   numstep is the number of steps
%   y0 is the initial condition

syms y(x)
cond = y(0) == y0;
A = dsolve(eqn, cond);
s1 = A;
S1 = s1;
for x = t0 : h : h*(numstep)
    subs(x);
    if x == t0    
        S1 = subs(s1,x);
    else 
      S1 = [subs(S1), subs(s1,vpa(x))];
    end
end
end

推荐阅读