Переписывал программу с R на Python, но почему-то значения y(x) на пайтоне слишком скачут в зависимости от N, при n=100 порядок ординат десятки, при 1000 порядок сотни, на R как было в пределах +- от 0 до 1, так и остается, получается какая - та неустойчивость решения на python. Кто-то может объяснить с чем мб связано дело?
Вот на R:
require(ggplot2)
ua = 0.8
ub = 1
A = function(x){
return( log(3 + sqrt(1 + x) ))
}
B = function(x){
return( exp(-x / 2) )
}
C = function(x){
return ( sqrt((3 - x) / 2) )
}
F = function(x){
return ( 2 - abs(1 - x) )
}
run = function(N){
h = 2.4/(N + 1)
an = function(x){
return( -(2 * A(x) + B(x) * h) )
}
bn = function(x){
return( 4 * A(x) + 2 * C(x) *h^2 )
}
cn = function(x){
return( -(2 * A(x) - B(x) * h) )
}
fn = function(x){
return( 2 * F(x) * h^2 )
}
x = seq(from = 0.1, to = 2.5, length.out = N+2)
u = vector("numeric", N+2)
a = sapply(x, an)
b = sapply(x, bn)
c = sapply(x, cn)
f = sapply(x, fn)
a[1] = 0
b[1] = 1
c[1] = 0
f[1] = 0.8
u[1] = 0.8
a[N+2] = 0
b[N+2] = 1
c[N+2] = 0
f[N+2] = 1
u[N+2] = 1
alpha = vector("numeric", N+2)
beta = vector("numeric", N+2)
alpha[1] = -c[1]/b[1]
beta[1] = (f[1])/b[1]
for(n in 2:(N + 1)){
alpha[n]= -c[n] / (a[n] * alpha[n-1] + b[n])
beta[n] = (f[n] - a[n] * beta[n-1]) / (a[n] * alpha[n-1] + b[n])
}
for(n in (N+1):2){
u[n] = alpha[n] * u[n+1] + beta[n]
}
df = data.frame(x = x, y = u)
g = ggplot(df, aes(x)) + geom_line(aes(y=y), col = "red")
print(g)
print(u[N+2])
}
run(100)
Вот на python:
import math as m
import matplotlib.pyplot as plt
def run(n):
h = 2.4 / (n+1)
a = list(range(n)); b = list(range(n)); c = list(range(n)); f = list(range(n))
# начальные условия
a[0] = 0; b[0] = 1; c[0] = 0 ; f[0] = 0.8
a[-1] = 0; b[-1] = 1; c[-1] = 0; f[-1] = 1
alpha = list(range(n)); beta = list(range(n))
alpha[0] = (-1)*(c[0])/(b[0])
beta[0] = (f[0])/(b[0])
x = list(range(n))
x[0] = 0.1
for i in range(1,n):
x[i] = x[0] + i*h
A = list(range(n)); B = list(range(n)); C = list(range(n)); F = list(range(n))
# функци A(x), B(x), C(x), f(x).
for i in range(n):
A[i] = m.log(3 + m.sqrt(1+x[i]))
B[i] = m.exp(-x[i]/2)
C[i] = m.sqrt((3-x[i])/2)
F[i] = 2 - m.fabs(1-x[i])
# коэф-ты.
for i in range(1,n):
a[i] = (-1) * (2 * A[i] + h * B[i])
b[i] = (4 * A[i] + 2 * pow(h, 2) * C[i])
c[i] = (-1) * (2 * A[i] - h * B[i])
f[i] = 2*pow(h, 2) * f[i]
alpha[i] = (-1)*(c[i])/(a[i]*alpha[i-1]+b[i])
beta[i] = (f[i] - a[i]*beta[i-1])/ (a[i]*alpha[i-1]+b[i])
# Решение u_(N+1)
y = list(range(n))
y[0] = 0.8
y[-1] = 1 # кладу последний n-й элемент равный 1
for i in reversed(range(1, n-1)):
y[i] = alpha[i]*y[i+1]+beta[i]
plt.plot(x, y, color='red', marker='o', linestyle='solid')
plt.show()
run(100)
