OCTAVE Solving Ordinary differential equation

Admin 1
Viewed 674

I have to solve the following differential equation with Octave

enter image description here enter image description here

Where

  • A is a matrix of (N-1)X(N-1) ,
  • f is defined as: f(x,t)=1000*sqrt(|1-t|) ,
  • μ=1,
  • u(x,0)=1000x(1-x)(1+(3/2)x^3),
  • the bar has lenght L=1 and u(0,t)=u(1,t)=0 for all t ,
  • h is 1/N (N number of partitions of the interval).

I have read how to use the command, but I do not understand how to use it in this specific case. How can I input the function du/dt. I'd be very grateful for your help.

differential-equations ode numerical-methods octave
1 Answers

I'm not an expert, and there's some information you've left out (e.g. the specific nature of x).

However, as an example of how to 'run' the lsode function, here's an example below, which makes the following assumption:

  • You mentioned a 'bar' of length 1, and N partitions/bins, so I'm assuming X is a vector representing the corresponding positions at the 'edges' of the bins, in which case you have N+1 edges. Since you are given the solutions for the 'outer' edges (i.e. x = 0 and x = 1), we only care for the 'inner edges', hence we only need to deal with vectors of N-1 elements (hence the size of matrix 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.' );

This may not be the same as your problem, but the above code runs as expected, so hopefully it will help you as a starting point.

One thing to note, the u vector passed into the lsode function needs to be in 'column' form. (the list of timepoints can be column or row, doesn't matter).

user_41831910
Related