程序如下
function
x=gauss(a,b)
%高斯求解方程组
%x=gauss(a,b)
n=length(a);
a=[a,b];
for
k=1:n-1
maxa=max(abs(a(k:n,k)));
if
maxa==0
return;
end
for
i=k:n
if
abs(a(i,k))==maxa
y=a(i,k:n+1);a(i,k:n+1)=a(k,k:n+1);a(k,k:n+1)=y;
break;
end
end
for
i=k+1:n
l(i,k)=a(i,k)/a(k,k);
a(i,k+1:n+1)=a(i,k+1:n+1)-l(i,k)a(k,k+1:n+1);
end
end
%回代
if
a(n,n)==0
return
end
x(n)=a(n,n+1)/a(n,n);
for
i=n-1:-1:1
x(i)=(a(i,n+1)-sum(a(i,i+1:n)x(i+1:n)))/a(i,i);
end
调用示例如下:
>>
a=[2,-1,3;4,2,5;1,2,0];
>>
b=[1;4;7];
>>
x=gauss(a,b)
x
=
9
-1
-6
昨天才回答过这个问题你可以再搜搜的
gauss消去法的分析。
其包括两个过程:
消去过程:把方程组系数矩阵a化为同解的上三角矩阵;
回代过程:按相反的顺序,从xn至x1逐个求解上三角方程组。
%高斯消去法的matlab程序
function
x=gauss(a,b);
%编写高斯消去法函数
%a表示方程组的系数矩阵,b表示方程组的值
%x表示最终的输出结果,即方程组的解
n=length(b);
%计算方程组的维数
%下面的程序在不断的消去,直到变成a变成上三角矩阵未知
for
k=1:n-1
for
i=k+1:n
a(i,k)=a(i,k)/a(k,k);
for
j=k+1:n
a(i,j)=a(i,j)-a(i,k)a(k,j);
end
b(i)=b(i)-a(i,k)b(k);
end
end
%表示高斯消去法的回带过程
x=zeros(n,1);
x(n)=b(n)/a(n,n);
for
k=n-1:-1:1
s=b(k);
for
j=k+1:n
s=s-a(k,j)x(j);
end
x(k)=s/a(k,k);
end
实例验证:
%调用编好的消去法函数
>>
a=[1,2,3;2,2,3;-1,-3,10];b=[0,3,2];gauss(a,b)
ans
=
30000
-15517
00345
>>
a=[1,2,3;2,2,3;-1,-3,10];b=[0,3,2];x=gauss(a,b)
x
=
30000
-15517
00345
>>
ax
%反代求解进行比较
ans
=
00000
30000
20000
高斯列主元消去法
function X=Gauss_pivot(A,B)
% 用Gauss列主主元消去法解线性方程组AX=B
%X是未知向量
n=length(B);
X=zeros(n,1);
c=zeros(1,n);
d1=0
for i=1:n-1
max=abs(A(i,i));
m=i;
for j=i+1:n
if max<abs(A(j,i))
max=abs(A(j,i));
m=j;
end
end
if(m~=i)
for k=i:n
c(k)=A(i,k);
A(i,k)=A(m,k);
A(m,k)=c(k);
end
d1=B(i);
B(i)=B(m);
B(m)=d1;
end
for k=i+1:n
for j=i+1:n
A(k,j)=A(k,j)-A(i,j)A(k,i)/A(i,i);
end
B(k)=B(k)-B(i)A(k,i)/A(i,i);
A(k,i)=0;
end
end
%回代求解
X(n)=B(n)/A(n,n);
for i=n-1:-1:1
sum=0;
for j=i+1:n
sum=sum+A(i,j)X(j);
end
X(i)=(B(i)-sum)/A(i,i);
End
用高斯消元法解线性方程组 的MATLAB程序
输入的量:系数矩阵 和常系数向量 ;
输出的量:系数矩阵 和增广矩阵 的秩RA,RB, 方程组中未知量的个数n和有关方程组解 及其解的信息
function [RA,RB,n,X]=gaus(A,b)
B=[A b]; n=length(b); RA=rank(A);
RB=rank(B);zhica=RB-RA;
if zhica>0,
disp('请注意:因为RA~=RB,所以此方程组无解')
return
end
if RA==RB
if RA==n
disp('请注意:因为RA=RB=n,所以此方程组有唯一解')
X=zeros(n,1); C=zeros(1,n+1);
for p= 1:n-1
for k=p+1:n
m= B(k,p)/ B(p,p); B(k,p:n+1)= B(k,p:n+1)-m B(p,p:n+1);
end
end
b=B(1:n,n+1);A=B(1:n,1:n); X(n)=b(n)/A(n,n);
for q=n-1:-1:1
X(q)=(b(q)-sum(A(q,q+1:n)X(q+1:n)))/A(q,q);
end
else
disp('请注意:因为RA=RB<n,所以此方程组有无穷多解')
end
end
1、下图是需要求解的线性方程组。
2、打开MATLAB,利用左除法(\)求解上述线性方程组。输入如下代码:close all; clear all; clc% MATLAB左除法(\)求解线性方程组,A = [1 2 3;-1 3 7;9 0 3];b = [1 4 7]';x = A\b。
3、保存和运行上述代码,利用左除法(\)得到线性方程组的解。
4、用求逆法(inv)求解线性方程组,输入如下代码:close all; clear all; clc,% MATLAB求逆法(inv)求解线性方程组,% A是线性方程组等号左边系数构成的矩阵。
5、保存和运行上述代码,利用求逆法(inv)得到线性方程组的解如下。
6、最后,可以看到左除法(\)和求逆法(inv)求得的线性方程组解是一样的。
大概就是个模拟,和平时解方程的过程、方法是一样的
program ttdd8;
var matrix:array[1100,1101] of double;
temp:array[1101] of double;
x:array[1100] of double;
i,j,k,n:integer;
m:double;
begin
readln(n);{一共有n个式子}
for i:=1 to n do
begin
for j:=1 to n+1 do
read(matrix[i,j]);{依次读入每个式子x1,x2…的系数及常数项}
readln;
end;
for i:=1 to n-1 do{需要消元n-1次}
begin
for j:=i to n do{寻找主元,即当前要消去元素系数最大的一个式子}
if matrix[j,i]>matrix[i,i] then
begin
temp:=matrix[i];
matrix[i]:=matrix[j];
matrix[j]:=temp;
end;
for j:=i+1 to n do
begin
m:=matrix[j,i]/matrix[i,i];
matrix[j,i]:=0;
for k:=i+1 to n+1 do
matrix[j,k]:=matrix[i,k]m-matrix[j,k];
end;
end;
x[n]:=matrix[n,n+1]/matrix[n,n];{回代过程}
for i:=n-1 downto 1 do
begin
m:=0;
for j:=i+1 to n do
m:=m+matrix[i,j]x[j];
x[i]:=(matrix[i,n+1]-m)/matrix[i,i];
end;
for i:=1 to n do
writeln(x[i]:0:2);
end
!高斯消去法
subroutine agaus(a,b,n,x,l,js)
dimension a(n,n),x(n),b(n),js(n)
double precision a,b,x,t
l=1 !逻辑变量
do k=1,n-1
d=00
do i=k,n
do j=k,n
if (abs(a(i,j))>d) then
d=abs(a(i,j))
js(k)=j
is=i
end if
end do
end do !把行绝对值最大的元素换到主元位置
if (d+10==10) then
l=0
else !最大元素为0无解
if(js(k)/=k) then
do i=1,n
t=a(i,k)
a(i,k)=a(i,js(k))
a(i,js(k))=t
end do !最大元素不在K行,K行
end if
if(is/=k) then
do j=k,n
t=a(k,j)
a(k,j)=a(is,j)
a(is,j)=t !交换到K列
end do
t=b(k)
b(k)=b(is)
b(is)=t
end if !最大元素在主对角线上
end if !消去
if (l==0) then
write(,100)
return
end if
do j=k+1,n
a(k,j)=a(k,j)/a(k,k)
end do
b(k)=b(k)/a(k,k) !求三角矩阵
do i=k+1,n
do j=k+1,n
a(i,j)=a(i,j)-a(i,k)a(k,j)
end do
b(i)=b(i)-a(i,k)b(k)
end do
end do
if (abs(a(n,n))+10==10) then
l=0
write(,100)
return
end if
x(n)=b(n)/a(n,n)
do i=n-1,1,-1
t=00
do j=i+1,n
t=t+a(i,j)x(j)
end do
x(i)=b(i)-t
end do
100 format(1x,'fail')
js(n)=n
do k=n,1,-1
if (js(k)/=k) then
t=x(k)
x(k)=x(js(k))
x(js(k))=t
end if
end do
return
end
program main
dimension a(4,4),b(4),x(4),js(4)
double precision a,b,x
real m1,m2,j
open(1,file="laiyitxt")
read(1,)m1,m2,j
close(1)
n=4
print,m1,m2,j
a(1,1)=m1cos(314159j/180)
a(1,2)=-m1
a(1,3)=-sin(314159j/180)
a(1,4)=0
a(2,1)=m1sin(314159j/180)
a(2,2)=0
a(2,3)=cos(314159j/180)
a(2,4)=0
a(3,1)=0
a(3,2)=m2
a(3,3)=-sin(314159j/180)
a(3,4)=0
a(4,1)=0
a(4,2)=0
a(4,3)=-cos(314159j/180)
a(4,4)=1
b(1)=0
b(2)=m198
b(3)=0
b(4)=m298
call agaus(a,b,n,x,l,js)
if (l/=0) then
write(,)"a1=",x(1),"a2=",x(2) ,"n1=",x(3),"n2=",x(4)
end if
end
!逆矩阵求解
SUBROUTINE qiuni(A,N,L,IS,JS)
DIMENSION A(N,N),IS(N),JS(N)
DOUBLE PRECISION A,T,D
L=1
DO K=1,N
D=00
DO I=K,N
DO J=K,N
IF(ABS(A(I,J))GTD) THEN !把最大的元素给D
D=ABS(A(I,J))
IS(K)=I
JS(K)=J
END IF
END DO
END DO
IF (D+10EQ10)THEN
L=0
WRITE(,200)
RETURN
END IF
200 FORMAT(1X,'ERRNOT INV')
DO J=1,N
T=A(K,J)
A(K,J)=A(IS(K),J)
A(IS(K),J)=T
END DO
DO I=1,N
T=A(I,K)
A(I,K)=A(I,JS(K))
A(I,JS(K))=T
END DO
A(K,K)=1/A(K,K)
DO J=1,N
IF(JNEK)THEN
A(K,J)=A(K,J)A(K,K)
END IF
END DO
DO I=1,N
IF(INEK)THEN
DO J=1,N
IF(JNEK)THEN
A(I,J)=A(I,J)-A(I,K)A(K,J)
END IF
END DO
END IF
END DO
DO I=1,N
IF(INEK)THEN
A(I,K)=-A(I,K)A(K,K)
END IF
END DO
END DO
DO K=N,1,-1
DO J=1,N
T=A(K,J)
A(K,J)=A(JS(K),J)
A(JS(K),J)=T
END DO
DO I=1,N
T=A(I,K)
A(I,K)=A(I,IS(K))
A(I,IS(K))=T
END DO
END DO
RETURN
END
SUBROUTINE BRMUL(A,B,N,C)
DIMENSION A(N,N),B(N),C(N)
DOUBLE PRECISION A,B,C
DO I=1,N
DO J=1,N
C(I)=00
DO L=1,N
C(I)=C(I)+A(I,L)B(L)
END DO
END DO
END DO
RETURN
END
program main
DIMENSION A(4,4),B(4,1),C(4,1),IS(4),JS(4)
DOUBLE PRECISION A,B,C
REAL M1,M2,JD
OPEN(1,FILE='LAIYITXT')
READ(1,) M1,M2,JD
PRINT,M1,M2,JD
CLOSE(1)
A(1,1)=M1COS(314JD/180)
A(1,2)=-M1
A(1,3)=-SIN(314JD/180)
A(1,4)=0
A(2,1)=M1SIN(314JD/180)
A(2,2)=0
A(2,3)=COS(314JD/180)
A(2,4)=0
A(3,1)=0
A(3,2)=M2
A(3,3)=-SIN(314JD/180)
A(3,4)=0
A(4,1)=0
A(4,2)=0
A(4,3)=-COS(314JD/180)
A(4,4)=1
B(1,1)=0
B(2,1)=M198
B(3,1)=0
B(4,1)=M298
CALL QIUNI(A,4,L,IS,JS)
CALL BRMUL(A,B,4,C)
WRITE(,) (C(I,1),I=1,4)
END
画图
USE MSFLIB
INTEGER status
TYPE(xycoord) xy
status=SETCOLORRGB(#FFFFFF)
status1=SETCOLORRGB(#0000FF)
OPEN(1,FILE="GTXT")
READ(1,) G1,G2,G3,G4
OPEN(2,FILE="NTXT")
READ(2,) N1,N2,N3,N4
CALL MOVETO(INT(20),INT(20),xy)
status=LINETO(INT(40),INT(G1))
status=LINETO(INT(80),INT(G2))
status=LINETO(INT(120),INT(G3))
status=LINETO(INT(160),INT(G4))
CALL SETLINESTYLE(#FF00)
CALL MOVETO(INT(20),INT(20),xy)
status1=LINETO(INT(40),INT(N1))
status1=LINETO(INT(80),INT(N2))
status1=LINETO(INT(120),INT(N3))
status1=LINETO(INT(160),INT(N4))
READ(,)
END
以上就是关于MATLAB用高斯消去法解非线性方程组的代码全部的内容,包括:MATLAB用高斯消去法解非线性方程组的代码、如何使用MATLAB,解答如下高斯消元法题目、找列主元高斯消去法来求解线性代数方程组解的matlab程序等相关内容解答,如果想了解更多相关内容,可以关注我们,你们的支持是我们更新的动力!
欢迎分享,转载请注明来源:内存溢出
评论列表(0条)