Main Content

本页的翻译已过时。点击此处可查看最新英文版本。

decic

ode15i 计算一致的初始条件

说明

示例

[y0_new,yp0_new] = decic(odefun,t0,y0,fixed_y0,yp0,fixed_yp0) 使用 y0yp0 作为完全隐函数 odefun 的初始条件的估计值,保留由 fixed_y0fixed_yp0 指定的分量作为固定分量,然后计算非固定分量的值。结果将获得一套完整的一致初始条件。新值 yo_newyp0_new 满足 odefun(t0,y0_new,yp0_new) = 0,适合用作 ode15i 的初始条件。

[y0_new,yp0_new] = decic(odefun,t0,y0,fixed_y0,yp0,fixed_yp0,options) 还使用 options 结构体 optionsAbsTolRelTol 指定值。可以使用 odeset 创建 options 结构体。

[y0_new,yp0_new,resnrm] = decic(___) 返回 odefun(t0,y0_new,yp0_new) 的范数作为 resnrm。如果范数看起来异常大,则使用 options 减小相对误差容限 RelTol,其默认值为 1e-3

示例

全部折叠

考虑隐式方程组

0=2y1-y20=y1+y2

这些方程很简单,很容易读取变量的一致初始条件。例如,如果固定 y1=1,则根据第二个方程 y2=-1,根据第一个方程 y1=-1/2。由于 y1y1y2 的这些值满足方程,因此它们是一致的。

通过使用 decic 计算方程的一致初始条件并固定值 y1=1,可以确认这些值。使用 y0 = [1 0]yp0 = [0 0] 的估计值,它们不满足方程,因此不一致。

odefun = @(t,y,yp) [2*yp(1)-y(2); y(1)+y(2)];
t0 = 0;
y0 = [1 0];
yp0 = [0 0];
[y0,yp0] = decic(odefun,t0,y0,[1 0],yp0,[])
y0 = 2×1

     1
    -1

yp0 = 2×1

   -0.5000
         0

计算一致的初始条件,并用 ode15i 求解隐式 ODE。

Weissinger 方程是

ty2(y)3-y3(y)2+t(t2+1)y-t2y=0

由于该方程采用一般形式 f(t,y,y)=0,您可以使用 ode15i 函数来求解隐式微分方程。

编写方程代码

要以适合 ode15i 的形式对方程进行编码,您需要编写一个函数,该函数使用 tyy 的输入并返回方程的残差值。可以使用函数 @weissinger 对此方程进行编码。查看函数文件。

type weissinger
function res = weissinger(t,y,yp)
%WEISSINGER  Evaluate the residual of the Weissinger implicit ODE
%
%   See also ODE15I.

%   Jacek Kierzenka and Lawrence F. Shampine
%   Copyright 1984-2014 The MathWorks, Inc.

res = t*y^2 * yp^3 - y^3 * yp^2 + t*(t^2 + 1)*yp - t^2 * y;

计算一致的初始条件

ode15i 求解器需要一致的初始条件,即提供给求解器的初始条件必须满足

f(t0,y,y)=0

由于可能会提供不一致的初始条件,并且 ode15i 不会检查一致性,因此建议您使用辅助函数 decic 来计算这些条件。decic 会固定一些指定的变量,并为不固定的变量计算一致的初始值。

在本例中,固定初始值 y(t0)=32,并让 decicy(t0)=0 的初始估计值开始,计算导数 y(t0) 的一致初始值。

t0 = 1;
y0 = sqrt(3/2);
yp0 = 0;
[y0,yp0] = decic(@weissinger,t0,y0,1,yp0,0)
y0 = 1.2247
yp0 = 0.8165

求解方程

使用 decic 加上 ode15i 返回的一致初始条件,在时间区间 [1 10] 上求解 ODE。

[t,y] = ode15i(@weissinger,[1 10],y0,yp0);

绘制结果

此 ODE 的精确解为

y(t)=t2+12

绘制 ode15i 计算的数值解 y 对解析解 ytrue 的图。

ytrue = sqrt(t.^2 + 0.5);
plot(t,y,'*',t,ytrue,'-o')
legend('ode15i', 'exact')

Figure contains an axes. The axes contains 2 objects of type line. These objects represent ode15i, exact.

输入参数

全部折叠

要求解的函数,指定为函数句柄,用于定义要求积分的函数。odefun 代表您要使用 ode15i 求解的隐式微分方程组。

对于标量 t 以及列向量 yyp 来说,函数 f = odefun(t,y,yp) 必须返回数据类型为 singledouble 的列向量 f,这些类型与 f(t,y,y') 相对应。odefun 必须同时接受 tyyp 这三个输入参数,即使其中某个参数未在函数中使用也是如此。

例如,要求解 y'y=0,请使用此函数。

function f = odefun(t,y,yp)
f = yp - y;

对于方程组,odefun 的输出为向量。每个方程变成解向量中的一个元素。例如,要求解

y'1y2=0y'2+1=0,

,请使用以下函数。

function dy = odefun(t,y,yp)
dy = zeros(2,1);
dy(1) = yp(1)-y(2);
dy(2) = yp(2)+1;

有关如何为函数 odefun 提供其他参数的信息,请参阅参数化函数

示例: @myFcn

数据类型: function_handle

初始时间,指定为标量。decic 使用初始时间计算满足 odefun(t0,y0_new,yp0_new) = 0 的一致初始条件。

数据类型: single | double

y 分量的初始估计值,指定为向量。y0 中的每个元素为 odefun 定义的方程组中的一个因变量 yn 指定一个初始条件。

数据类型: single | double

要保持固定的 y 分量,指定为向量 1 和 0,或者指定为 []

  • 如果 y0(i) 的估计值中不允许进行任何更改,则设置 fixed_y0(i) = 1

  • 如果可以更改任何条目,则设置 fixed_y0 = []

您最多只能固定 length(yp0) 个分量。根据具体问题的不同,可能无法始终固定 y0yp0 的某些分量。最好不要固定不必要的多余分量。

y' 分量的初始估计值,指定为向量。yp0 中的每个元素为 odefun 定义的方程组中的一个微分因变量 y'n 指定一个初始条件。

数据类型: single | double

要保持固定的 y' 分量,指定为向量 1 和 0,或者指定为 []

  • 如果 yp0(i) 的估计值中不允许进行任何更改,则设置 fixed_yp0(i) = 1

  • 如果可以更改任何条目,则设置 fixed_yp0 = []

您最多只能固定 length(yp0) 个分量。根据具体问题的不同,可能无法始终固定 y0yp0 的某些分量。最好不要固定不必要的多余分量。

Options 结构体,指定为结构体数组。使用 odeset 函数创建或修改 options 结构体。可与 decic 函数一起使用的相关选项包括 RelTolAbsTol,它们控制用于计算初始条件的误差阈值。

示例: options = odeset('RelTol',1e-5)

数据类型: struct

输出参数

全部折叠

y0 的一致初始条件,以向量形式返回。如果 resnrm 的值较小,则 yo_newyp0_new 满足 odefun(t0,y0_new,yp0_new) = 0,适合用作 ode15i 的初始条件。

yp0 的一致初始条件,以向量形式返回。如果 resnrm 的值较小,则 yo_newyp0_new 满足 odefun(t0,y0_new,yp0_new) = 0,适合用作 ode15i 的初始条件。

残差范数,以向量形式返回。resnrmodefun(t0,y0_new,yp0_new) 的范数。

  • 如果 resnrm 的值较小,则表示 decic 已成功计算满足 odefun(t0,y0_new,yp0_new) = 0 的一致初始条件。

  • 如果 resnrm 的值很大,请使用 options 输入尝试调整误差阈值 RelTolAbsTol

提示

  • 在使用 ode15i 求解之前,ihb1daeiburgersode 示例文件使用 decic 计算一致初始条件。可以键入 edit ihb1daeedit iburgersode 查看代码。

  • 还可以使用 decic 计算通过 ode15sode23t 求解的 DAE 的一致初始条件。为此,请按以下步骤操作。

    1. 采用完全隐式的形式 f(t,y,y') = 0 重写方程组。

    2. 调用 decic 计算方程的一致初始条件。

    3. 指定 y0_new 作为调用求解器的初始条件,并指定 yp_new 作为 odesetInitialSlope 选项的值。

另请参阅

| |

在 R2006a 之前推出