作者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