Chebyshev Differentiation Matrix to solve ODE

75 views (last 30 days)
Abel
Abel on 16 Oct 2012
Answered: Francis Oketch on 30 Jul 2020
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!

Answers (3)

Star Strider
Star Strider on 18 Oct 2012
Edited: Star Strider on 19 Oct 2012
It took longer than I thought it would for me to learn something about spectral methods and particularly Chebyshev differentiation matrices. Fortunately, the relevant chapters of Spectral Methods in MATLAB are available online.
The problem you are having seems to be in the way you define your f variable:
f = 2*x.*ones(size(x));
in the context of your original differential equation:
dy/dx = 2*x*y
Specifically, y is not a vector of ones. (I am not sure how you would define it in that situation.) When I defined f as:
f = 2*x.*exp(x.^2 - 1);
and then calculated:
u = lsqr(D1,f);
I got a value for u that differed from the analytic integrated differential equation by a constant of 0.6773 but otherwise matched it with acceptable precision. (I haven't accounted for that particular constant of integration.)
I also calculated u using lsqr because using the mldivide appropriate for dense matrices threw a condition number warning when applied to this problem.

Qiqi Wang
Qiqi Wang on 11 May 2017
Edited: Qiqi Wang on 11 May 2017
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;
  2 Comments
Lukgaf
Lukgaf on 19 Mar 2018
I used the cod but not calable. kindly help
Walter Roberson
Walter Roberson on 19 Mar 2018
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.

Sign in to comment.


Francis Oketch
Francis Oketch on 30 Jul 2020
% Solve: y'=2*x*exp(x^2-1); y(-1)=1
clc; clear all;
N=20;
% compute the chebyshev differentiation matrix and x-grid
[D,x] = cheb(N);
% compute the derivatives at x (i.e. at the chebyshev grid points)
A=D;
R = 2*x.*exp(x.^2 - 1);
% boundary condition (correct way)
A(end,:) = 0;
A(end,end) = 1;
R(end) = 1;
% solve
u = A\R;
plot(x,u,'--rs',x,exp(x.^2-1))
legend('Chebyshev','Exact')

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!