matlabpde

How do I change the poisson equation to solve for an alternative answer in Matlab


I'm trying to redo the Poisson equation so that ∂²Φ /∂x² = 2x² - 0.5x + exp(x). Every time I try and input the right-hand side of the equation it comes up with a syntax error, any help would very much be appreciated.

% Solving of 1D Poisson equation
% using finite differences

% Clearing memory and figures
clear all
clf

% Define 1D numerical model
xsize=100; % Model size in horizontal direction, m
Nx=101; % Number of grid points
dx=xsize/(Nx-1); % Gridstep, m
x=0:dx:xsize; % Coordinates of grid points, m

% Defining global matrixes
L=zeros(Nx,Nx); % Koefficients in the left part
R=zeros(Nx,1); % Right parts of equations


% Composing global matrixes
% by going through all grid points
for j=1:1:Nx
    % Discriminating between BC-points and internal points
    if(j==1 || j==Nx)
        % BC-points
        L(j,j)=1; % Left part
        R(j)=0; % Right part
    else
        % Composing Poisson eq. d2T/dx2=1
        %
        % ---T(j-1)------T(j)------T(j+1)---- STENCIL
        %
        % 1/dx2*T(j-1)+(-2/dx2)*T(j)+1/dx2*T(j+1)=1
        %
        % Left part of j-equation
        L(j,j-1)=1/dx^2; % T(j-1)
        L(j,j)=-2/dx^2; % T(j)
        L(j,j+1)=1/dx^2; % T(j+1)
        % Right part
        R(j)=2x;
    end
end

% Solve the matrix
S=L\R;

% Reload solutions
T=zeros(1,Nx); % Create array for temperature, K
for j=1:1:Nx
    T(j)=S(j);
end

% Visualize results
figure(1); clf
plot(x,T,'o r') % Plotting results

Again, I appreciate any help in this matter


Solution

  • The problem is clear: "input the right-hand side of the equation it comes up with a syntax error", to correct the error you should properly define the right-hand side of the equation:

    R(j)=2*x(j)^2-0.5*x(j) + exp(x(j));
    

    Since it depends on x, you should point to each element j when calculating the coefficients in the for loop on the left and right hand sides of the equation.

    Analytical solution of the PDE is: enter image description here

    Changing this line of code and plotting the discretized and analytical solutions produces

    % Visualize results
    figure(1); clf
    plot(x,T,'o r');
    hold on
    plot(x , 0.166667*x.^4-0.0833333*x.^3-268811714181613573321934397213431784538112*x+exp(x)-1)
    

    enter image description here

    which seems correct. Differences are due to discretization itself.