Issue
I am trying to use piecewise
to integrate function that change at certain domain.
The result of the 3D wireframe plot (with spb) is different than when using manual way of computing integral one by one. The correct plot is the one that is using manual way.
If you are wondering, why there is for loop from 1
to 10
, it is to replace sum from 1 to infinity, since it can represent sum of 1 to infinity (the Fourier Series in Heat Equation). It converges anyway in the end.
This is the plot (the correct plot is on the right):
So, where is my fault in this piecewise MWE:
# https://stackoverflow.com/questions/77259323/can-sympy-compute-definite-integral-by-cases/77259532#77259532
# https://stackoverflow.com/questions/77258199/how-to-replace-the-variable-from-sympy-computation-so-it-can-be-plotted-with-mat/77258592#77258592
import numpy as np
import sympy as sm
from sympy import *
from spb import *
x = sm.symbols("x")
t = sm.symbols("t")
n = sm.symbols("n", integer=True)
L = 20
D = 0.475
f = (S(2)/L)*sin(n*np.pi*x/20)*Piecewise((x, (0 <= x) & (x <= 10)), (20-x, (10 < x) & (x <= 20)))
print('The function u(x,0) : ')
print('')
sm.pretty_print(f)
print('')
print('')
print(piecewise_fold(f))
fint = integrate(f, (x, 0, 20))
g = fint*exp(-(n**2)*(np.pi**2)*D*t/400).nsimplify()
print('')
print('')
sm.pretty_print(fint)
print('')
print('')
s3 = 0
for c in range(10):
s3 += g.subs({n:c})
print('')
print('The function u(x,t) : ')
print('')
sm.pretty_print(s3)
plot3d(
s3, (x, 0, 20), (t, 0, 10), {"alpha": 0}, # hide the surface
wireframe=True, wf_n1=20, wf_n2=10,
wf_rendering_kw={"color": "tab:blue"}, # optional step to customize the wireframe lines
backend=MB, zlabel="$u(x,t)$", title="One Dimensional Heat Equation"
)
The code with the correct plot, but still integrating with manual way, not using piecewise
:
# https://stackoverflow.com/questions/77258199/how-to-replace-the-variable-from-sympy-computation-so-it-can-be-plotted-with-mat/77258592#77258592
import numpy as np
import sympy as sm
from sympy import *
from spb import *
x = sm.symbols("x")
t = sm.symbols("t")
n = sm.symbols("n", integer=True)
L = 20
f1 = (2/L)*x*sin(n*np.pi*x/20)
f2 = (2/L)*(20-x)*sin(n*np.pi*x/20)
fint1 = sm.integrate(f1,(x,0,10))
fint2 = sm.integrate(f2,(x,10,20))
D = 0.475
g = (fint1+fint2)*sin(n*np.pi*x/20)*exp(-(n**2)*(np.pi**2)*D*t/400).nsimplify()
s = 0
for c in range(10):
s += g.subs({n:c})
print(s)
print('')
print('The function u(x,t) : ')
print('')
sm.pretty_print(s)
print('')
print('')
plot3d(
s, (x, 0, 20), (t, 0, 10), {"alpha": 0}, # hide the surface
wireframe=True, wf_n1=20, wf_n2=10,
wf_rendering_kw={"color": "tab:blue"}, # optional step to customize the wireframe lines
backend=MB, zlabel="$u(x,t)$", title="One Dimensional Heat Equation"
)
Solution
The discrepancy comes from this:
g = fint*exp(-(n**2)*(np.pi**2)*D*t/400).nsimplify()
and this:
g = (fint1+fint2)*sin(n*np.pi*x/20)*exp(-(n**2)*(np.pi**2)*D*t/400).nsimplify()
Note that in the second one you have multiplied by sin(n*np.pi*x/20)
, whereas the first one is missing that multiplier.
On a side note, when using SymPy always use SymPy objects. So, in your code you want to replace np.pi
with pi
. It also makes code easier to read.
Answered By - Davide_sd
0 comments:
Post a Comment
Note: Only a member of this blog may post a comment.