MHB Plotting an infinite domain PDE

AI Thread Summary
The discussion focuses on plotting an infinite domain partial differential equation (PDE) with a solution expressed as an integral. Users are seeking efficient methods to visualize this in Matlab, Mathematica, or Python, noting challenges with computation time, particularly in Mathematica. A provided Mathematica code integrates over limited bounds, while a Python implementation struggles with accuracy despite handling larger integration ranges. Suggestions include optimizing the integration process by factoring out constants and recognizing that the inner integral may not converge, potentially utilizing the Riemann-Lebesgue Lemma for further analysis. The conversation emphasizes the need for more accurate and efficient plotting techniques for complex PDE solutions.
Dustinsfl
Messages
2,217
Reaction score
5
Is anyone familiar with plotting an infinite domain PDE where the solution is an integral.

Take the solution
\[
T(x,t) = \frac{100}{\pi}\int_0^{\infty}\int_{-\infty}^{\infty} \frac{\sinh(u(10-y)}{\sinh(10u)} \cos(u(\xi-x))d\xi du
\]
How could I plot this in Matlab, Mathematica, or Python?

As a note, Mathematica is slow because I have code that does it, but when I increase u, the time to complete goes up exponentially. I am talking hours.

Everything I tried in Matlab failed.
 
Physics news on Phys.org
I'm curious: can you post your Mathematica code?
 
Ackbach said:
I'm curious: can you post your Mathematica code?

Code:
f[x_, y_, u_, \[Xi]_] =  Sinh[u*(5 - y)]/Sinh[u*5]*100/\[Pi]*Cos[u*(\[Xi] - x)];
 
g[x_?NumericQ, y_?NumericQ] := NIntegrate[f[x, y, u, \[Xi]], {u, 1, 10}, {\[Xi], -5, 5}];

DensityPlot[g[x, y], {x, -5, 5}, {y, 0, 5}, PlotPoints -> 5]

This code only integrate u from 1 to 10 and \(\xi\) from -5 to 5. However, I would like to increase the integration to have a more accurate plot.

Here is Python code that takes a long time but can handle a bigger integration bounds in the same time Mathematica can but the code isn't producing an accurately plot either so there must be something wrong as well.

Code:
import numpy as np
import pylab as pl 
from scipy.integrate import dblquad  

b = 50.0  
x = np.linspace(-10, 10, 1000) 
y = np.linspace(0, 10, 1000)  def f(xx, u, xi, yj):     
    return ((np.exp(u * (b - yj)) - np.exp(-u * (b - yj))) / 
              (np.exp(u * b) - np.exp(-u * b)) * np.cos(u * (xx - xi)))  T = np.zeros([len(x), len(y)])

 
for i, xi in enumerate(x):     
    for j, yj in enumerate(y):         
        T[i, j] += dblquad(             
                   f, -10, 10, lambda x: 0.1, lambda x: 10, args=(xi, yj))[0]fig = pl.figure()
ax = fig.add_subplot(111)
pax = ax.pcolormesh(y, x, T)
fig.colorbar(pax)
pl.xlim((-X, X))
pl.ylim((0, b))
pl.grid(True)

pl.show()
 
You might save some computational time if you pull everything out of the inner integral that you can:
$$T(x,t)= \frac{100}{ \pi} \int_{0}^{ \infty} \frac{ \sinh(u(10-y))}{ \sinh(10u)}
\underbrace{ \int_{- \infty}^{ \infty} \cos(u( \xi-x)) \, d \xi}_{ \text{define this as }f(u,x)} \, du.$$
This way, the ratio on sinh's isn't evaluated for every inner loop. I should point out, however, that the inner integral, which is $f(u,x)$, doesn't converge. The Riemann-Lebesgue Lemma might come in handy here, depending on where you want to look at things.
 

Similar threads

Replies
4
Views
3K
Replies
1
Views
4K
Replies
3
Views
3K
Replies
1
Views
5K
Back
Top