时间戳

首页 > 解决方案 > 求解具有分段系数的 ODE

matlab - 求解具有分段系数的 ODE

问题描述

我必须在整个区间上绘制[0,500]二阶 ODE 的解y''(x)+f(x)y(x)+y(x)=0，初始条件y(0)=0.1为 , y'(0)=0.1, 其中f : [0,500] -> [0,1],

在此处输入图像描述

我使用了这个功能

function dwdt=equat(t,w)
f = @(x) zeros(size(x));
for i= 2:500
    Bounds_1 = [i - (1/i^2),i];
    Bounds_2 = [i,i + (1/i^2)];
    f1 = @(x) (i^2.*x - i^3 + 1).*(Bounds_1(1) < x & x <= Bounds_1(2));
    f2 = @(x) (-i^2.*x + i^3 + 1).*(Bounds_2(1) < x & x <= Bounds_2(2));
    f = @(x) f(x) + f1(x) + f2(x); 
end
dwdt=zeros(2,1);
dwdt(1)=w(2);
dwdt(2)=-f*w(2)-w(1);
end

用命令

tspan = [0 500];
z0=[0.1 0.1];
[t,z] = ode45(@(t,z) equat(t,z), tspan, z0);
figure
plot(t,z(:,1),'r');

我收到几个错误

Unary operator '-' is not supported for operand of type 'function_handle'.

Error in equat (line 12)
dwdt(2)=-f*w(2)-w(1);

Error in @(t,z)equat(t,z)

Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)
  odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

我怎样才能修复它，以获得解决方案的绘图y(x)？

LE：使用 Lutz 建议的函数绘制解决方案：在此处输入图像描述

标签： matlabode

解决方案

所以你在i半径的位置有“帽子”或三角形尖峰1/i^2。由于它们不重叠i>=2，因此只有一个是活动的，非零。因此你可以计算

i = round(x); i = max(2,min(500,i));
fx = max(0,1-i*i*abs(x-i));

如果不通过匿名函数进行大量扭曲，这是无法实现的，但是您可以将其直接包含到 ODE 系统函数中或将其声明为普通函数。

function fx = f(x)
  i = max(2,min(500,round(x)));
  fx = max(0,1-i*i*abs(x-i));
end%function

function dwdt=equat(x,w)
  dwdt = zeros(2,1);
  dwdt(1) = w(2);
  dwdt(2) = -f(x)*w(2)-w(1);
end%function

您可能会观察到，积分器可能会跳过后面的非常细的尖峰。更严格地设置误差容限，设置最大步长（这可能非常不经济）或将分段线性函数拆分为线性的部分，如给定的那样，并用于odextend将各个部分拼接在一起。

推荐阅读