Chebyshev Differentiation Matrix to solve ODE
조회 수: 36 (최근 30일)
이전 댓글 표시
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!
댓글 수: 0
채택된 답변
Shrinivas Chimmalgi
2024년 4월 4일
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
댓글 수: 1
John D'Errico
2024년 4월 4일
편집: John D'Errico
2024년 4월 4일
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.
추가 답변 (1개)
Qiqi Wang
2017년 5월 11일
편집: Qiqi Wang
2017년 5월 11일
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
Walter Roberson
2018년 3월 19일
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.
참고 항목
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!