MATLAB用高斯消去法解非线性方程组的代码

MATLAB用高斯消去法解非线性方程组的代码,第1张

程序如下

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程序等相关内容解答,如果想了解更多相关内容,可以关注我们,你们的支持是我们更新的动力!

欢迎分享,转载请注明来源:内存溢出

原文地址: http://outofmemory.cn/zz/9984292.html

(0)
打赏 微信扫一扫 微信扫一扫 支付宝扫一扫 支付宝扫一扫
上一篇 2023-05-04
下一篇 2023-05-04

发表评论

登录后才能评论

评论列表(0条)

保存