Welcome to our community

Be a part of something great, join today!

Plotting an infinite domain PDE

dwsmith

Well-known member
Feb 1, 2012
1,673
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.
 

Ackbach

Indicium Physicus
Staff member
Jan 26, 2012
4,192
I'm curious: can you post your Mathematica code?
 

dwsmith

Well-known member
Feb 1, 2012
1,673
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()
 

Ackbach

Indicium Physicus
Staff member
Jan 26, 2012
4,192
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.