Chebyshev Differentiation Matrix to solve ODE

Trying to solve dy/dx = 2xy, y(1) = 1, using Chebyshev differentiation matrix

The exact solution is y = exp(x.^2 -1);

Heres what I have:

% compute the chebyshev differentiation matrix and x-grid
[D1,x] = cheb(N);

% boundary condition (I don't know if this is correct?)
D1 = D1(1:end-1,1:end-1); x = x(1:end-1);

% compute the derivatives at x (i.e. at the chebyshev grid points)
f = 2*x.*ones(size(x));

% solve 
u = D1\f;

% set the boundary condition
u = [u;1];

Where cheb.m is from Trefethen (spectral methods in matlab)

function [D,x] = cheb(N)

    % check the base case
    if N == 0; D = 0; x = 1; return; end

    % create the Chebyshev grid
    x = cos(pi*(0:N)/N)';

    c = [2; ones(N-1,1);2].*(-1).^(0:N)';
    X = repmat(x,1,N+1);
    dX = X-X';
    D = (c*(1./c)')./(dX+(eye(N+1)));
    D = D - diag(sum(D'));

This solution (u = D1\f) does not match the exact solution at all.

I think what I have is close ... Any help would be awesome. Thanks in advance!

Answer 1

[Shrinivas Chimmalgi]

Although this question is more than a decade old, there still seems to be some interest in this. I was suprised to see that all the answers use the solution y = exp(x.^2 -1) which is what we are supposed to find by solving the ODE dy/dx = 2xy, y(1) = 1 numerically.

clear
clc
% Solve: dy/dx = 2xy, y(x), y(1) = 1
N = 20;
% compute the chebyshev differentiation matrix and x-grid
[D,x] = cheb(N); %dom=[-1,1]
% Domain is important as we will be solving for y only inside this domain.

% The ODE problem can now be written as the following linear problem with
% one constraint. This constraint arises due to the boundary condition.
% D*y_est == 2*x.*y_est
% y(x=1) = 1.

y_est = nan(N+1,1);
y_est(end) = 1; % boundary condition

% Note that we have N unknowns (values of y at x(1:N)) but the system above has N+1 equations.
% To solve this we utilize the boundary condition (constraint) to reduce
% the system to the form q = Ap with size(A) = [N,N] (N equations).
% Solving this using mldivide (\ operator) gives us
y_est(1:N) = (D(1:N,1:N)-2*diag(x(1:N)))\(-D(1:N,N+1)*1);

% Remark that we didn't use y = exp(x.^2 -1). In the figure below we can
% see that the solution to the ODE we calculated matches the exact solution
% y = exp(x.^2 -1) well at the x-grid points. 
% Plotting
plot(x,y_est,'s',x,exp(x.^2-1))
xlabel("x")
ylabel("y, y_{est}")
legend('Chebyshev','Exact')

%%
% Checking convergence. Here we can see how the numerically computed
% solution improves as we increase number of x-grid points N.
figure
for N = 2:30
% compute the chebyshev differentiation matrix and x-grid
[D,x] = cheb(N); %dom=[-1,1]
% Domain is important as we will be solving for y only inside this domain.

% The ODE problem can now be written as the following linear problem with
% one constraint. This constraint arises due to the boundary condition
% y(x=1) = 1.

% D*y_est == 2*x.*y_est
y_est = nan(N+1,1);
y_est(end) = 1; % boundary condition

% Note that we have N unknowns (values of y at x(1:N)) but the system above has N+1 equations.
% To solve this we utilize the boundary condition (constraint) to reduce
% the system to the form q = Ap with size(A) = [N,N].
% Solving this using mldivide (\ operator) gives us
y_est(1:N) = (D(1:N,1:N)-2*diag(x(1:N)))\(-D(1:N,N+1)*1);

semilogy(N,norm(y_est-exp(x.^2-1)),'bo')
xlabel('N')
ylabel("||y-y_{est}||")
hold on
drawnow
end

function [D,x] = cheb(N)
% check the base case
if N == 0; D = 0; x = 1; return; end
% create the Chebyshev grid
x = cos(pi*(0:N)/N)';
c = [2; ones(N-1,1);2].*(-1).^(0:N)';
X = repmat(x,1,N+1);
dX = X-X';
D = (c*(1./c)')./(dX+(eye(N+1)));
D = D - diag(sum(D'));
end

[John D'Errico]

I had to laugh myself when I was reading through some of the other answers, since I was observing the same thing. People seem to be trying to solve an ode using the known exact solution. The good news is, I could accept your answer.

Answer 2

[Qiqi Wang]

The original post was right that the boundary condition was not set correct. Here's a right way to set the boundary condition:

% compute the chebyshev differentiation matrix and x-grid
[D,x] = cheb(N);
% compute the derivatives at x (i.e. at the chebyshev grid points)
f = 2*x.*exp(x.^2 - 1);
% boundary condition (correct way)
D(end,:) = 0;
D(end,end) = 1;
f(end) = 1;
% solve 
u = D\f;

[Lukgaf]

I used the cod but not calable. kindly help

[Walter Roberson]

What value are you passing? What error message are you receiving?

You opened a Question on this topic, so the discussion should probably be there.