简介
如果你想知道宏伟的计划...阅读简介。如果没有,请跳到我的问题。
我有一个微分方程和线性代数课程的项目,我必须使用计算机代数系统来解决具有常系数的一阶线性常微分方程和具有常系数的非线性常微分方程。我必须通过分析和数字来证明这是完成的。我的计划是有 2 个函数,这些函数利用dsolve
函数作为分析部分,以及一个名为ODE1
的函数,这是我通过 Matlab 的本教程学到的。该函数被定义为使用欧拉方法近似微分方程的解。这些函数都将使用相同的输入参数,因此每个输入可以定义一次,并且这些函数都将了解要使用的输入(可能将函数嵌套在一个调用函数下(。目标是评估用于查找数值近似值的相同区间的解析解,并在表格和图形中比较结果。我得到了数值解决方案,给我一个行向量形式的"表格",并绘制该行向量。我开始遇到分析解决方案的问题...
我的问题
为了求解一阶线性常微分方程,我生成了这个函数
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,当使用 0:0.4 步在区间内使用 5 个步骤进行评估时,列表应该是
1
1.105...
1.2214...
1.3498...
1.4918...
经过一些四舍五入。
所以这里的问题是我的分析解决方案中有一个额外的
1
条目。我相信它与for
循环中S1 = [subs(S1); vpa(subs(s1))]
的subs(S1)
部分有关,但我对如何解决这个问题感到困惑。
我有点理解为什么我需要使用subs
函数,因为我正在使用符号变量来使用在其答案中输出符号变量的dsolve
函数。此外,为了使for
循环迭代和更改,每次都必须将符号变量替换为t
的实际值。我确实尝试在for
循环之前移动vpa(subs(s1))
,但这只是在向量中返回了 5 次相同的值。我也尝试不使用subs(S1)
,它给了我
exp(t)
1.0
1.1051709...
1.2214027...
1.3498588...
所以我肯定这是代码的这一部分。
旁注:我知道分析方法输出列向量,链接的视频中显示的ODE1
也是如此。为了让 Matlab 将其绘制为一条线,我转置了列向量以制作行向量,并在解部分固定后对解析解执行相同的操作。
通过更改for
循环的内部结构,我使它工作。我的最终函数代码是这样的: 函数 [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