Нашел ошибку. Исправил:
http://narod.ru/disk/start/02.dl3b-narod.yandex.ru/38456229001/h6011cd5ec3e741387bb28e9ccf15430d/PDE_exam.pdfНапишу здесь код на питоне
(Оффтоп)
from scipy import *
from scipy.integrate import quad,dblquad
from pylab import *
a=sqrt(1.1)*0.01
x=0.5
b=2500
pi=3.14
def U1(t):
res, error = dblquad(lambda tau, xi:
10.0/sqrt(t-tau)*( exp(-( (xi-x)**2.0/a**2.0/4.0/(t-tau) ))
- exp(-( (xi+x)**2.0/a**2.0/4.0/(t-tau) )) ),
0, Inf, lambda xi: 0, lambda xi: t)
res *= 1.0/2.0/a/sqrt(pi)
return res
def U2(t):
res, error = quad(lambda tau:
b*exp(-t)/(t-tau)*exp(-x**2.0/4/a**2.0/(t-tau)), 0, t)
res *= a/sqrt(pi)
return res
t1=0.1
t2=300.0
n=30
results1 = linspace(t1, t2, n)
results2 = linspace(t1, t2, n)
results = linspace(t1, t2, n)
i=0
for item in linspace(t1, t2, n):
results1[i] = U1(item)
results2[i] = U2(item)
results[i] = results1[i] + results2[i]
print "t = ", item, " U(0.5, t) = ", results[i]
i+=1
print ">>>>> U(0.5, 200) = ", U1(200.0)+U2(200.0)
figure()
scatter(linspace(t1, t2, n), results)
show()
И результаты
(Оффтоп)
t = 0.1 U(0.5, t) = 1.0342252742e-18
t = 10.4413793103 U(0.5, t) = 104.44026986
t = 20.7827586207 U(0.5, t) = 207.880286113
t = 31.124137931 U(0.5, t) = 311.320302358
t = 41.4655172414 U(0.5, t) = 414.760314387
t = 51.8068965517 U(0.5, t) = 518.20022473
t = 62.1482758621 U(0.5, t) = 621.640557217
t = 72.4896551724 U(0.5, t) = 725.076407655
t = 82.8310344828 U(0.5, t) = 828.502672156
t = 93.1724137931 U(0.5, t) = 931.908471449
t = 103.513793103 U(0.5, t) = 1035.27797882
t = 113.855172414 U(0.5, t) = 1138.58841209
t = 124.196551724 U(0.5, t) = 1241.82017601
t = 134.537931034 U(0.5, t) = 1344.94469192
t = 144.879310345 U(0.5, t) = 1447.9341081
t = 155.220689655 U(0.5, t) = 1550.76016612
t = 165.562068966 U(0.5, t) = 1653.39497386
t = 175.903448276 U(0.5, t) = 1755.81159026
t = 186.244827586 U(0.5, t) = 1857.98444033
t = 196.586206897 U(0.5, t) = 1959.88958848
t = 206.927586207 U(0.5, t) = 2061.50489891
t = 217.268965517 U(0.5, t) = 2162.81010925
t = 227.610344828 U(0.5, t) = 2263.78684068
t = 237.951724138 U(0.5, t) = 2364.42330202
t = 248.293103448 U(0.5, t) = 2464.69547692
t = 258.634482759 U(0.5, t) = 2564.59477738
t = 268.975862069 U(0.5, t) = 2664.1097831
t = 279.317241379 U(0.5, t) = 2763.23047908
t = 289.65862069 U(0.5, t) = 2861.94824058
t = 300.0 U(0.5, t) = 2960.25567512
>>>>> U(0.5, 200) = 1993.46690129
Результат Mathcad:
U(x=0.5, t=200) = 50460
Графики похожие, если не считать маткадовский более вогнутым. Маткадовский ответ в 25 раз больше питоновского. Где ошибка? Кто может в своих программах посчитать?