OCTAVE 求解常微分方程

Question

我必须用 Octave 解决以下微分方程


在哪里

A是 (N-1)X(N-1) 的矩阵,

f定义为:f(x,t)=1000*sqrt(|1-t|) ,

μ=1 ,

u(x,0)=1000x(1-x)(1+(3/2)x^3) ,

酒吧有长度L=1和 u(0,t)=u(1,t)=0所有 t ,

h是 1/N ( N 区间的分区数)。

我已经阅读了如何使用该命令，但我不明白在这种特定情况下如何使用它。如何输入函数 du/dt。我将非常感谢您的帮助。

Answer 1

我不是专家，您遗漏了一些信息(例如 x 的特定性质)。
但是，作为如何“运行” lsode 函数的示例，下面是一个示例，它做出以下假设:

您提到了一个长度为 1 的“条”和 N 个分区/分箱，因此我假设 X 是一个向量，表示分箱“边缘”处的相应位置，在这种情况下，您有 N+1 条边。由于您得到了“外”边(即 x = 0 和 x = 1)的解，我们只关心“内边”，因此我们只需要处理 N-1 个元素的向量(因此大小矩阵 A)。

  N         = 4;   % number of bins
  mu        = 1;
  h         = 1 / N;
  A_matrix  = rand( N - 1, N - 1 );
  t_rvector = [0 : 20].';

  x_cvector = linspace( 0, 1, N + 1 ).';   % N+1 edges, including outer ones
  x_cvector = x_cvector(2 : end - 1);      % discard outer edges, keep inner only.


% Create f(t) as defined in the question, but vectorised w.r.t. x
% (i.e. outputting an N-1 element vector that can be added to a corresponding u vector)

  f_function = @(t) 1000 .* sqrt( abs( 1 - t ) ) .* ones( size( x_cvector ) );


% Create a u0 vector (of N-1 elements, corresponding to inner edges on the bar)
  u0_cvector = 1000 .* x_cvector .* (1 - x_cvector) .* (1 + (3 / 2) .* x_cvector .^ 3);

% Create a 'wrapper' function of the form expected by `lsode`
  ode_function = @(u_cvector, t) -mu ./ (h .^ 2) .* A_matrix * u_cvector + f_function(t);

% Run `lsode`
  [u_matrix, ISTATE, MSG] = lsode( ode_function, u0_cvector, t_rvector )

% Pad u_matrix with 0s on left and right, corresponding to x=0 and x=1 values
  TimePoints = length( t_rvector );
  ZeroColumn = zeros( TimePoints, 1 );
  u_matrix   = [ ZeroColumn, u_matrix, ZeroColumn ];

% Plot the curves (when plotting a matrix, each column is treated as a separate plot), therefore if you want to plot along the 'components' of u (corresponding to positions in $
  plot ( u_matrix.' );

这可能与您的问题不同，但上面的代码按预期运行，因此希望它可以帮助您作为起点。
需要注意的一件事， u矢量传递到 lsode功能需要采用“列”形式。 (时间点列表可以是列或行，无关紧要)。