Giải hệ phương trình

[pththao]

Code đó tương ứng với phương trình này:

Mình có một hệ hai phương trình như sau:
x = K(a-x)(b-x-y)
y =K'(b-y)(b-x-y)
Trong đó, x và y là ẩn, còn lại K, K', a, b là hằng số.
Bạn nào có kinh nghiệm có thể cho biết là về nguyên tắc thì có giải được hệ trên không?

Thông qua thuật toán này:

Về nguyên tắc phương trình trên tương đương với phương trình bậc ba nên giải được tổng quát theo công thức Cardano, hoặc dùng mathematica là nhanh nhất:
Đưa hệ về dạng:
t = b - x - y = x/(K(a-x)) = y/(K'(b-y)),
trong đó ta đặt biến mới t. Từ hai hệ thức cuối giải được x và y theo t,
x = aKt/(aKt+1), y = bK't/(bK't+1).
Thế x và y vào phương trình đầu,
t = b - aKt/(Kt+1) - bK't/(K't+1),
hay
a + t = b/(K't+1) + a/(Kt+1).
Phương trình trên sau khi qui đồng mẫu số sẽ là bậc ba. Cũng có thể có thủ thuật nào giải nhanh phương trình đó nhưng mình chưa nghĩ ra

Nhìn qua thấy phương trình bạn vừa post khác với phương trình nêu ban đầu?
Thống nhất pt nào đi, mai mình đến chỉnh thử code xem.

 

[Ho Huu Tho]

Mình cứ nghĩ hai hệ giống nhau, tương ứng với code mà Thảo cung cấp. Thống nhất là dùng hệ mà mình vừa post xong nhé:
x = K1(a-x)(a-x-y)
y = K2(b-y)(a-x-y)

Mình nghĩ hai hệ giống nhau vì:
Theo gợi ý của Thảo, mình đặt ẩn phụ là t = a-x-y, ta có:
x=aK1t/(K1t+1)
y=bK2t/(K2t+1)
Mình thấy hai cái phương trình này giống hệt cái trong bộ code của Thảo:

#Solving:
t0 = itersolve(t0 ,a ,b , K1, K2, delta)
x = a*K1*t0/(K1*t0+1)
y = b*K2*t0/(K2*t0+1)

 

[pththao]

OK, mai mình xem lại. Thêm một ý nhỏ, phương trình trên chưa mô tả sự bắt mồi của P' với A', (P'+A' = P'A'). Mình nghĩ giải hệ ba ẩn đó cũng không khó khăn gì lắm (nếu quá trình lặp hội tụ, không thì chỉ có computer là khổ tí thôi.) Thọ muốn mô tả chính xác hơn thì viết thử hệ ba phương trình đó lúc rảnh mình code luôn xem sao.

 

[Ho Huu Tho]

Cảm ơn Thảo đã gợi ý, đúng là bỏ qua phương trình còn lại sẽ không chính xác.
Hệ ba phương trình phản ứng như sau:
A + A' = AA'; hệ số cân bằng là K1
p + A' = PA'; hệ số cân bằng là K2
P' + A = P'A; hệ số cân bằng là K3

Nồng độ ban đầu của A = A' là a; của P là b; của P' là c.
Lúc cân bằng có nồng độ của AA' là x; của PA' là y; của P'A là z.

-> nồng độ lúc cân bằng của A là a-x-z, của A' là a-x-y, của P là b-y, của P' là c-z.

Ta có hệ phương trình (hệ 3) cho hệ lúc cân bằng như sau:
x = K1(a-x-z)(a-x-y)
y = K2(b-y)(a-x-y)
z = K3(c-z)(a-x-z)

 

[pththao]

Mình hỏi thêm, cỡ giá trị thông thường của K1, K2, K3, điều này liên quan đến tính hội tụ của quá trình lặp...

Nhân tiện, pt trên có khác, một phương trình là -a một phương trình là -b.

 

[pththao]

Mình cứ nghĩ hai hệ giống nhau, tương ứng với code mà Thảo cung cấp. Thống nhất là dùng hệ mà mình vừa post xong nhé:
x = K1(a-x)(a-x-y)
y = K2(b-y)(a-x-y)

Mình nghĩ hai hệ giống nhau vì:
Theo gợi ý của Thảo, mình đặt ẩn phụ là t = a-x-y, ta có:
x=aK1t/(K1t+1)
y=bK2t/(K2t+1)
Mình thấy hai cái phương trình này giống hệt cái trong bộ code của Thảo:

#Solving:
t0 = itersolve(t0 ,a ,b , K1, K2, delta)
x = a*K1*t0/(K1*t0+1)
y = b*K2*t0/(K2*t0+1)

Mình chỉnh viết lại một đoạn đơn giản hơn, dùng tỷ lệ x/a, y/b làm đầu ra cho biến. Rất tiếc điều kiện là đầu vào không được thay đổi quá xa để quá trình là hội tụ Với giá trị Thọ cho vào quá trình của V[2] thổi ra vô hạn nên hết đường cứu chữa. Với các đặt biến khéo léo thế nào đó sẽ xử lý được, nhưng mình chưa thấy phương hướng ở đâu.

#HISTORY:
#Febuary 24, 2011: Created.
#Febuary 26, 2011: Modified to vector itterative form.
rm(list = ls())
#Function iteratively solve an eq.
#V0 = estimation of solution.
#a,K1,b,K2: parameters (rescaled).
#delta: error threshold.
#return values: vector of 2 components, (x/a,y/b).
itersolve = function(V0, a, b, K1, K2, delta = 0.0001)
{
N = 100000
V1 = V0
V2 = V0
BAD = TRUE

for (k in 1:N)
{
V2[1] = K1*(1-V1[1])*(a - a*V1[1] - b* V1[2])
V2[2] = K2*(1-V1[2])*(a - a*V1[1] - b* V1[2])
V = V1-V2
s = sqrt(V %*% V)
if (s > delta)
{
V1 = V2
}
else
{
BAD = FALSE
break
}
}
if (BAD)
{
stop('Solution is badly approximated, function exited!')
}
return(V2)
}
#This can be called from the cmd.
itersolve(V0 = c(0,0), a = 1, b = 1, K1 = 0.01, K2 = 0.02, delta = 0.0001)

 

[pththao]

Cảm ơn Thảo đã gợi ý, đúng là bỏ qua phương trình còn lại sẽ không chính xác.
Hệ ba phương trình phản ứng như sau:
A + A' = AA'; hệ số cân bằng là K1
p + A' = PA'; hệ số cân bằng là K2
P' + A = P'A; hệ số cân bằng là K3

Nồng độ ban đầu của A = A' là a; của P là b; của P' là c.
Lúc cân bằng có nồng độ của AA' là x; của PA' là y; của P'A là z.

-> nồng độ lúc cân bằng của A là a-x-z, của A' là a-x-y, của P là b-y, của P' là c-z.

Ta có hệ phương trình (hệ 3) cho hệ lúc cân bằng như sau:
x = K1(a-x-z)(a-x-y)
y = K2(b-y)(a-x-y)
z = K3(c-z)(a-x-z)

Code tương ứng, (không biết có nhìn nhầm chỗ nào không.) Cũng như trên không cho vào được các giá trị quá xa của K1, K2, K3 (yêu cầu nhỏ hơn 1), a0, b0, c0 không quá lớn

#HISTORY:
#Febuary 24, 2011: Created.
#Febuary 26, 2011: Modified to vector itterative form.
#Febuary 26, 2011: Modified to system of 3 eqs.
rm(list = ls())
#Function iteratively solve an eq.
#V0 = estimation of solution.
#a,K1,b,K2: parameters (rescaled).
#delta: error threshold.
#return values: vector of 2 components, (x/a,y/b).
itersolve = function(V0, a0, b0, c0, K1, K2, K3, delta = 0.0001)
{
N = 100000
V1 = V0
V2 = V0
BAD = TRUE

for (k in 1:N)
{
V2[1] = K1/a0*(a0 - a0*V1[1] - c0*V1[3])*(a0 - a0*V1[1] - b0*V1[2])
V2[2] = K2/b0*(b0 - b0*V1[2])*(a0 - a0*V1[1] - b0*V1[2])
V2[3] = K3/c0*(c0 - c0*V1[3])*(a0 - a0*V1[1] - c0*V1[3])
V = V1-V2
s = sqrt(V %*% V)
if (s > delta)
{
V1 = V2
}
else
{
BAD = FALSE
break
}
}
if (BAD)
{
stop('Solution is badly approximated, function exited!')
}
return(V2)
}
#This can be called from the cmd.
itersolve(V0 = c(0,0,0), a0 = 1, b0 = 1, c0 = 1, K1 = 0.01, K2 = 0.02, K3 = 0.001, delta = 0.0000001)

 

[pththao]

Sau khi source('file.R'), bạn có thể giọi liên tiếp bằng lệnh từ màn hình (không phái source tiếp nữa). Vi dụ:

>source('file.R')
>itersolve(V0 = c(0,0,0), a0 = 1, b0 = 1, c0 = 1, K1 = 0.01, K2 = 0.02, K3 = 0.001, delta = 0.0000001)

Phương pháp xấp xỷ liên tiếp như sau:
1. Mọi hệ phương trình đưa được về dạng X = f(X); trong đó X có thể là vector, ở trên là (x,y) và (x,y,z) cho hai trường hợp.
2. Thực hiện dãy:
X_{0} = xấp xỉ ban đầu.
X_{n+1} = f(X_{n})
3. Nếu dãy hội tụ đến một giá trị nào đó thì giá trị đó là nghiệm của phương trình, nếu không tạm thời kết luận xấp xỷ ban đầu kém. Ở trên khi a, b tăng quá quá trình không hội tụ. Với một biến đổi hợp lý, chẳng hạn X = tY ta có thể biến quá trình không hội tụ thành hội tụ, nhưng mình chưa nhìn ra chỗ ấy.
Điều kiện hội tụ phụ thuộc bản thân hàm f và giá trị ban đầu.

 

[Ho Huu Tho]

Ta có hệ phương trình (hệ 3) cho hệ lúc cân bằng như sau:
x = K1(a-x-z)(a-x-y)
y = K2(b-y)(a-x-y)
z = K3(c-z)(a-x-z)

Khoảng giá trị của a, b, c lần lượt là: 1e-10, 1e-7, 1e-7
Khoảng giá trị của K1, K2, K3 lần lượt là: 1e20, 1e7, 1e7

Nếu dãy không hội tụ thì code của Thảo sẽ báo đầu ra thế nào?

Nếu đặt ẩn phụ: x' = xK1; y' = yK2; z' = zK3, bằng cách biến đổi như sau liệu có được ích lợi gì không Thảo:

Nhân hai vế của từng phương trình của hệ 3 lần lượt với K1, K2, K3, ta có:
xK1 = (aK1-xK1-zK1)(aK1-xK1-yK1)
yK2 = (bK2-yK2)(aK2-xK2-yK2)
zK3 = (cK3-zK3)(aK3-xK3-zK3)

Gọi: x' = xK1; y' = yK2; z' = zK3, thay vào ta được hệ sau (hệ 4):
x' = (aK1-x'-K1*z'/K3)(aK1-x'-K1*y'/K2)
y' = (bK2-y')(aK2-K2*x'/K1-y')
z' = (cK3-z')(aK3-K3*x'/K1-z')

Nếu thay hệ phương trình này vào hệ phương trình vào bộ code của Thảo thì còn phải thay đổi gì khác nữa không? Mình copy đoạn code của Thảo và thay các phương trình mới bằng chữ màu đỏ như ở dưới:

#HISTORY:
#Febuary 24, 2011: Created.
#Febuary 26, 2011: Modified to vector itterative form.
#Febuary 26, 2011: Modified to system of 3 eqs.
rm(list = ls())
#Function iteratively solve an eq.
#V0 = estimation of solution.
#a,K1,b,K2: parameters (rescaled).
#delta: error threshold.
#return values: vector of 2 components, (x/a,y/b).
itersolve = function(V0, a0, b0, c0, K1, K2, K3, delta = 0.0001)
{
N = 100000
V1 = V0
V2 = V0
BAD = TRUE

for (k in 1:N)
{
V2[1] = (K1*a0 - V1[1] - K1*V1[3]/K3)*(K1*a0 - V1[1] - K1*V1[2]/K2)
V2[2] = (K2*b0 - V1[2])*(K2*a0 - K2*V1[1]/K1 - V1[2])
V2[3] = (K3*c0 - V1[3])*(K3*a0 - K3*V1[1]/K1 - V1[3])
V = V1-V2
s = sqrt(V %*% V)
if (s > delta)
{
V1 = V2
}
else
{
BAD = FALSE
break
}
}
if (BAD)
{
stop('Solution is badly approximated, function exited!')
}
return(V2)
}
#This can be called from the cmd.
itersolve(V0 = c(0,0,0), a0 = 1, b0 = 1, c0 = 1, K1 = 0.01, K2 = 0.02, K3 = 0.001, delta = 0.0000001)

 

[pththao]

Ta có hệ phương trình (hệ 3) cho hệ lúc cân bằng như sau:
x = K1(a-x-z)(a-x-y)
y = K2(b-y)(a-x-y)
z = K3(c-z)(a-x-z)

Khoảng giá trị của a, b, c lần lượt là: 1e-10, 1e-7, 1e-7
Khoảng giá trị của K1, K2, K3 lần lượt là: 1e20, 1e7, 1e7

Nếu dãy không hội tụ thì code của Thảo sẽ báo đầu ra thế nào?

Khi không hội tụ có hai trường hợp báo lỗi: 1. Lỗi mình báo, stop('Solution is badly approximated, function exited!'), 2. Lỗi máy báo do quá trình nhanh quá đẫn đến V = infinity sớm, thường là lỗi không định nghĩa phép toán nào đó (chẳng hạn so sánh.)

Nếu đặt ẩn phụ: x' = xK1; y' = yK2; z' = zK3, bằng cách biến đổi như sau liệu có được ích lợi gì không Thảo:

Nhân hai vế của từng phương trình của hệ 3 lần lượt với K1, K2, K3, ta có:
xK1 = (aK1-xK1-zK1)(aK1-xK1-yK1)
yK2 = (bK2-yK2)(aK2-xK2-yK2)
zK3 = (cK3-zK3)(aK3-xK3-zK3)

Gọi: x' = xK1; y' = yK2; z' = zK3, thay vào ta được hệ sau (hệ 4):
x' = (aK1-x'-K1*z'/K3)(aK1-x'-K1*y'/K2)
y' = (bK2-y')(aK2-K2*x'/K1-y')
z' = (cK3-z')(aK3-K3*x'/K1-z')

Nếu thay hệ phương trình này vào hệ phương trình vào bộ code của Thảo thì còn phải thay đổi gì khác nữa không? Mình copy đoạn code của Thảo và thay các phương trình mới bằng chữ màu đỏ như ở dưới:

#HISTORY:
#Febuary 24, 2011: Created.
#Febuary 26, 2011: Modified to vector itterative form.
#Febuary 26, 2011: Modified to system of 3 eqs.
rm(list = ls())
#Function iteratively solve an eq.
#V0 = estimation of solution.
#a,K1,b,K2: parameters (rescaled).
#delta: error threshold.
#return values: vector of 2 components, (x/a,y/b).
itersolve = function(V0, a0, b0, c0, K1, K2, K3, delta = 0.0001)
{
N = 100000
V1 = V0
V2 = V0
BAD = TRUE

for (k in 1:N)
{
V2[1] = (K1*a0 - V1[1] - K1*V1[3]/K3)*(K1*a0 - V1[1] - K1*V1[2]/K2)
V2[2] = (K2*b0 - V1[2])*(K2*a0 - K2*V1[1]/K1 - V1[2])
V2[3] = (K3*c0 - V1[3])*(K3*a0 - K3*V1[1]/K1 - V1[3])
V = V1-V2
s = sqrt(V %*% V)
if (s > delta)
{
V1 = V2
}
else
{
BAD = FALSE
break
}
}
if (BAD)
{
stop('Solution is badly approximated, function exited!')
}
return(V2)
}
#This can be called from the cmd.
itersolve(V0 = c(0,0,0), a0 = 1, b0 = 1, c0 = 1, K1 = 0.01, K2 = 0.02, K3 = 0.001, delta = 0.0000001)

Theo mình thấy thay đổi như thế đúng rồi, nguyên tắc chỉ đơn giản vậy thôi. Việc đổi biến thường có động lực từ bản chất vật lý và hóa học của bài toán: ví dụ có a = 10^10 mol chất tham gia phản ứng, con số đó rõ ràng lớn vô nghĩa, vì vậy ta định nghĩa tỷ lệ chất tham gia phản ứng x. Khi đó x thuộc khoảng [0,1] và trong thực tế là mong đợi x có giá trị thu nhận được cỡ 0.1 hoặc tương tự thế. Thọ đổi biến như thế trong trường hợp này cũng có thể có lợi, cái đó phải thử thôi, nhưng mình không thấy động lực vật lý để đổi biến như vậy