I am trying to calculate and plot subsonic coefficients of pressure symbolically, but it is returning a piecewise function that I can't easily plot, what am I doing wrong?

I am trying to calculate and plot subsonic coefficients of pressure symbolically, but it is returning a piecewise function that I can't easily plot, am I doing something?

clear all; clc; close all

syms x ksi

%N = 20;
gamma = 1.4;
limits = [0 1];
M = 1.2;
%alpha = deg2rad([-5:1:20]);
%X = linspace(0,1,N);

%Geometry
c = 1;
tau = 1/10; %thickness ratio
ZksiU = 2*tau*c*(ksi/c-(ksi/c)^2);
ZksiL = -0.5*tau*c*(ksi/c-(ksi/c)^2);

fun = diff(ZksiU,ksi,1)/(x-ksi);
Cpi = -(2/pi)*int(fun,[limits(1) limits(2)])    %incompressible pressure coefficient
CpSub = -(2/M^2*(gamma+1))*((1-M^2)-((1-M^2)^3/2+(3/4)*M^2*(gamma+1)*Cpi)^(2/3))    %Subsonic pressure coefficient 

returns:

fun = ((2*ksi)/5 - 1/5)/(ksi - x)

Cpi = piecewise(ksi in Dom::Interval([0], [1]), (5734161139222659*((2*ksi)/5 - 1/5)*int(1/(x - ksi), x, 0, 1))/9007199254740992, ~in(ksi, 'real') | ksi < 0 | 1 < ksi, -(5734161139222659*((2*ksi)/5 - 1/5)*(log(-ksi) - log(1 - ksi)))/9007199254740992)

CpSub = piecewise(ksi in Dom::Interval([0], [1]), (10*((464467052277035379*int(-1/(ksi - x), x, 0, 1)*((2*ksi)/5 - 1/5))/281474976710656000 - 1331/31250)^(2/3))/3 + 22/15, ~in(ksi, 'real') | ksi < 0 | 1 < ksi, (10*(- (464467052277035379*((2*ksi)/5 - 1/5)*(log(-ksi) - log(1 - ksi)))/281474976710656000 - 1331/31250)^(2/3))/3 + 22/15)

Answer 1

[Walter Roberson]

syms x ksi

%N = 20;
gamma = 1.4;
limits = [0 1];
M = 1.2;
%alpha = deg2rad([-5:1:20]);
%X = linspace(0,1,N);

%Geometry
c = 1;
tau = 1/10; %thickness ratio
ZksiU = 2*tau*c*(ksi/c-(ksi/c)^2)

ZksiU = 

ZksiL = -0.5*tau*c*(ksi/c-(ksi/c)^2);

fun = diff(ZksiU,ksi,1)/(x-ksi)

fun = 

Cpi = -(2/pi)*int(fun,[limits(1) limits(2)])    %incompressible pressure coefficient

Cpi = 

CpSub = -(2/M^2*(gamma+1))*((1-M^2)-((1-M^2)^3/2+(3/4)*M^2*(gamma+1)*Cpi)^(2/3))    %Subsonic pressure coefficient 

CpSub = 

KSI = linspace(-1, 2, 1000);
CP = double(subs(CpSub, ksi, KSI));
plot(KSI, real(CP), 'k-.', KSI, imag(CP), 'r-+')

vpa(subs(CpSub, 0.5))

ans = 

NaN

Why is it getting NaN for a fair bit of the region? Well, you are asking to integrate 1/(x-ksi) with ksi in the range 0 to 1, and x in the range 0 to 1. Except for the case of ksi = 0 exactly or ksi = 1 exactly, you are trying to integrate through a fundamental discontinuity, equivalent to integrating through 1/X with X being 0 at some point in the range.

Hence, the integral is only defined for ksi outside the 0 to 1 range.