作者swallowyann (歸心似燕)
看板MATLAB
標題[問題]非線性方程組，求聯立方程所有跟
時間Tue Jul 2 18:54:37 2013
各位大大晚安！
最近為了求解非線性方程組而傷透腦筋，本人有參考一本書【應用數值分析】，其中有提
到運用高斯消去法及牛頓迭代法可以解出，但是我運用其code卻無法跑出解，請各位大大
幫忙解題了！
(1)gauss99.m (2)A_nonlinear_newton.m (3)func1.m %func1是放我要解的聯立方程組
(1)
function [x]=gauss99(A,b)
%本方法為Guass消去法
%A=input('輸入矩陣A=');
%b=input('輸入矩陣b=');
[m,n]=size(A);
exa=[A b];
%fprintf('原方程式係數為:\n')
gg=0;
for i=1:n
for j=1:n+1
% fprintf('%8.3f ',exa(i,j));
end
% fprintf(';\n');
end
%scaling
for i=1:n
max=exa(i,1);
for j=2:n
if abs(max)<abs(exa(i,j))
max=exa(i,j);
end
end
for jj=1:n+1
if abs(max)<0.00001
gg=1;exa(i,jj)=0;
else
exa(i,jj)=exa(i,jj)/max;
end
end
end
%fprintf('scaling後方程式係數為:\n')
for k=1:n-1
%以下為轉軸 Pivoting
jjj=0;
max=exa(k,k);
for kk=k+1:n
if abs(max)<abs(exa(kk,k))
max=exa(kk,k);
jjj=kk;
end
end
if jjj~=0
if k~=jjj
for i=1:n+1
temp=exa(k,i);
temp1=exa(jjj,i);
exa(k,i)=temp1;
exa(jjj,i)=temp;
end
end
end
for i=k+1:n
if abs(exa(k,k))<0.00001
gg=1;pp=0;
else
pp=exa(i,k)/exa(k,k);
end
for j=k:n+1
exa(i,j)=exa(i,j)-pp*exa(k,j);
end
end
end
%fprintf('轉軸及消去後方程式係數為:\n')
for i=1:n
for j=1:n+1
end
end
%以下為往回代入法 backward substitute
if abs(exa(n,n))<0.0001
gg=1;
x(n)=999;
else
x(n)=exa(n,n+1)/exa(n,n);
end
for i=n-1:-1:1
rr=0;
for k=i+1:n
rr=rr+exa(i,k)*x(k);
end
if gg==1
x(n)=999;
x(i)=999;
else
x(i)=(exa(i,n+1)-rr)/exa(i,i);
end
end
(2)
function A_nonlinear_newton(func1,gauss99)
%解非線性聯立方程式
n=input('輸入方程式數目n=');
ii=0;ik=0;
xb=input('範圍下界xb=');
xn=input('範圍上界xn=');
tol=input('輸入允許誤差tol=');
max=input('輸入最大計算次數max=');
nn=(xn-xb)/0.5;
for i1=1:nn+1
xk(1)=xb+(i1-1).*0.5;
for i2=1:nn+1
xk(2)=xb+(i2-1).*0.5;
for i3=1:nn+1
xk(3)=xb+(i3-1).*0.5;
ii=0;
for i=1:n
x(i)=xk(i);
end
while (1)
for i=1:n
for j=1:n
y(j)=x(j);
if i==j
y(j)=x(j)+0.001;
end
end
z=feval(func1,x);
zz=feval(func1,y);
for k=1:n
d(k,i)=(zz(k)-z(k))/0.001;
end
end
if x(1)==999
break;
end
A=[d];
b=-[z];
for i6=1:n
for j6=1:n
if abs(A(i6,j6))>100000
break;
end
end
end
for i6=1:n
if abs(b(i6))>100000
break;
end
end
bb=b';
[xx]=feval(gauss99,A,bb);
pp=1;
for i=1:n
if abs(xx(i))>999
pp=0;
xx(i)=999;
break;
end
if abs(xx(i))>tol
pp=0;
end
end
if pp==1
ik=ik+1;
for i=1:n
xy(i,ik)=x(i)+xx(i);
exx(i,ik)=xx(i);
g(i,ik)=xk(i);
end
break;
end
ii=ii+1;
if ii>max
break;
end
for i=1:n
x(i)=x(i)+xx(i);
end
end
end
end
end
ka=1;
for i=1:n
fprintf('第%d組答案:x(%d)=%f;誤差x(%d)=%f',ka,i,xy(i,1),i,exx(i,1));
fprintf(',初始值座標x(%d)=%f\n',i,g(i,1));
ans(i,ka)=xy(i,1);
err(i,ka)=xy(i,1);
gk(i,ka)=g(i,1);
end
for j=2:ik
pi=0;
%for k=1:ka
k=0;
while (k<ka)
k=k+1;
% pi=0;
for i=1:n
if abs(xy(i,j)-ans(i,k))<0.01
pi=pi+1;
end
end
if pi==n
pi=1;
break;
else
pi=0;
end
end
if pi==0
ka=ka+1;
for ip=1:n
ans(ip,ka)=xy(ip,j);err(ip,ka)=exx(ip,j);gk(ip,ka)=g(ip,j);
fprintf('第%d組答案:x(%d)=%f,誤差x(%d)=%f',ka,ip,ans(ip,ka),ip,err(ip,ka));
fprintf(',初始值座標x(%d)=%f;\n',ip,gk(ip,ka))
end
end
end
(3)
function y = func1(x)
%定義非線性方程式_FULL_PARAMETER
P = 1;
Q1 = 4;
Q2 = 3;
Q3 = 2;
N0 = 3;
N1 = 2;
N2 = 1;
N3 = 1;
T = 42;
SINR = 5;
%L21 = 3.6055;
%L22 = 4.5826;
L21 = 3;
L22 = 4;
L23 = 5;
c1 = 4*((2*Q3)+Q2)*Q3*N1*T;
c2 = 2*SINR*Q1*N0*(2*Q3+Q2+Q1)*Q3;
c3 = (-2)*P*(2*Q3+Q2+Q1)*Q3*N1*T;
c4 = (-2)*(P^2)*((2*Q3+Q2)^3)*((N1*T)^2);
c5 = (-SINR)*Q1*N0*(2*Q3+Q2+Q1)*P*(2*Q3+Q2);
c6 = (P^2)*(2*Q3+Q2+Q1)*(2*Q3+Q2)*N1*T;
c7 = (-P)*2*Q3*N2*T;
c8 = (-SINR)*P*Q2*N0*(2*Q3+Q2);
c9 = (N0/(N1*T))*2*Q3*(2*Q3+Q2+Q1)*N2*T;
c10 = -(N0/(N1*T))*P*(2*Q3+Q2+Q1)*(2*Q3+Q2)*N2*T;
c11 = (P^2)*N2*T;
c12 = (-SINR)*((2*Q3)/(2*Q3+Q2+Q1))*(N0/N3*T)*((2*L21)+(2*L22)+(2*L23));
c13 = N0/(N1*T);
c14 = (P/(2*Q3+Q2+Q1))-(N0/(N2*T));
y(1) = (c1*x(1)) + (c2*x(2)) - (c3*x(1)*x(3)) - (c5*x(2)*x(3)) +
(c6*x(1)*x(2)*x(3));
y(2) = -(c7*x(1)*x(2)) + (c8*x(1)*x(3)) + (c9*x(2)) - (c10*x(2)*x(3)) +
(c11*x(1)*x(2)*x(3));
y(3) = -(c12*x(1)) - (c13*x(2)*x(3)) + (c14*x(1)*x(2)*x(3));
------------------------------------------------------------------------------
以上為3個.檔，執行後顯示以下訊息
>> A_nonlinear_newton('func1','gauss99')
輸入方程式數目n=3
範圍下界xb=0
範圍上界xn=10
輸入允頂~差tol=0.001
輸入最大計算次數max=20
??? Undefined function or variable "xy".
Error in ==> A_nonlinear_newton at 86
fprintf('第%d組答案:x(%d)=%f;誤差x(%d)=%f',ka,i,xy(i,1),i,exx(i,1));
------------------------------------------------------------------------------
我有想過可能是因為區間設定的問題，但試過好多區間依然這麼顯示，而因為我要求解的
值為正值，所以現在不知道是code的問題還是有其它問題存在，麻煩大大們解題了！
--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 60.244.175.119
1F：→ swallowyann:忘了補充最後要求的是x(1)、x(2)、x(3)，三個數值 07/02 18:57
2F：→ swallowyann:x(3)理應較大，x(1)最小！ 07/02 18:58
3F：→ jengfu:樓主大大＠＠ 真的有人會幫你看這麼長的文章嗎? 07/22 20:57