MATLAB计算方法算法代码
- ❤️❤️1. LU分解
- 💙函数文件.m
- 💚测试脚本.m
- 💜测试结果
- ❤️❤️2. Gauss消元
- 💙函数文件.m
- 💚测试脚本.m
- 💜测试结果
- ❤️❤️3. 平方根法
- 💙函数文件.m
- 💚测试脚本.m
- 💜测试结果
- ❤️❤️4. Jacobi迭代
- 💙函数文件.m
- 💚测试脚本.m
- 💜测试结果
- ❤️❤️5. Gauss-Seidel迭代
- 💙函数文件.m
- 💚测试脚本.m
- 💜测试结果
- ❤️❤️6. 正交多项式曲线拟合
- 💙函数文件.m
- 💚测试脚本.m
- 💜测试结果
- ❤️❤️7. 龙贝格积分
- 💙函数文件.m
- 💚测试脚本.m
- 💜测试结果
💞💞💞💞💞💞💞💞
💦💦💦💦💦💦💦💦
💫💫💫💫💫💫💫💫
❤️❤️1. LU分解
💙函数文件.m
function x= LU(A,b)
format short
%此函数用于利用LU分解法解方程组
r = rank(A);
[m,n]=size(A);
% [L,U,P]=lu(A);
% b = P*b;
L=eye(n);
for i=1:n-1
for j=i+1:n
L(j,i)=A(j,i)/A(i,i);
A(j,:)=A(j,:)-(A(j,i)/A(i,i))*A(i,:);
end
end
U=A;
L
U
%解Ly=b
y(1) = b(1);
if m>1
for i=2:m
y(i)=b(i)-L(i,1:i-1)*y(1:i-1)';
end
end
y = y';
%解Ux=y得方程组的解
x0(r)=y(r)/U(r,r);
if r>1
for i=r-1:-1:1
x0(i)=(y(i)-U(i,i+1:r)*x0(i+1:r)')/U(i,i);
end
end
x0 = x0';
if flag==1
x=x0;
return;
else
format rat;
Z = null(A,'r');
[mz,nz]=size(Z);
x0(r+1:n)=0;
for i=1:nz
t=sym(char([107 48+i]));
k(i)=t;
end
x=x0;
for i=1:nz
x = x+k(i)*Z(:,i);
end
end
💚测试脚本.m
clear;
clc;
format short
fprintf('示例：\n');
a = [2 2 3;4 7 7;-2 4 5]
b = [3 1 -7]'
result = LU(a,b)
fprintf('======================================================================================================\n\n');
💜测试结果
示例：
a =
2 2 3
4 7 7
-2 4 5
b =
3
1
-7
L =
1 0 0
2 1 0
-1 2 1
U =
2 2 3
0 3 1
0 0 6
result =
2
-2
1
======================================================================================================
>>
❤️❤️2. Gauss消元
💙函数文件.m
function [result] = gauss(a1,b)
% gauss 此处显示有关此函数的摘要
% a1-系数矩阵,b-右端常数列
[m,n] = size(a1); % 获取矩阵的阶数
a = [a1,b]; % 定义增广矩阵
for i = 1 : n-1 % i 每个主元素所在行
for j = i + 1 : n
if a(i,i) ~= 0
m = -a(j,i)/a(i,i); % 定义初等行变换所需的倍数
a(j,:) = a(j,:) + m * a(i,:); % 初等行变换，以得到约化主元下方的0元
else
disp('算法失败')
break
end
end
fprintf('********************第%.0f步消元结束如下：********************\n',i)
disp(a)
% pause
end
% 回代
b = a(:,n+1)
result(n) = b(n) / a(n,n);
fprintf('方程组的解result(%.0f) = %.5f：\n',n,result(n))
for k = (n-1):-1:1
s = 0;
for t = k + 1 : n
s = s + a(k,t) * result(t);
end
result(k) = (b(k)-s)/a(k,k);
fprintf('方程组的解result(%.0f) = %.5f\n',k,result(k));
end
💚测试脚本.m
clear;
clc;
fprintf('示例：\n');
a = [1 0 2 4;4 3 0 1;5 2 1 6;7 4 3 2]
b = [5 6 7 8]'
result = gauss(a,b)
fprintf('======================================================================================================\n\n');
💜测试结果
示例：
a =
1 0 2 4
4 3 0 1
5 2 1 6
7 4 3 2
b =
5
6
7
8
********************第1步消元结束如下：********************
1 0 2 4 5
0 3 -8 -15 -14
0 2 -9 -14 -18
0 4 -11 -26 -27
********************第2步消元结束如下：********************
1 0 2 4 5
0 3 -8 -15 -14
0 0 -11/3 -4 -26/3
0 0 -1/3 -6 -25/3
********************第3步消元结束如下：********************
1 0 2 4 5
0 3 -8 -15 -14
0 0 -11/3 -4 -26/3
0 0 0 -62/11 -83/11
b =
5
-14
-26/3
-83/11
方程组的解result(4) = 1.33871：
方程组的解result(3) = 0.90323
方程组的解result(2) = 4.43548
方程组的解result(1) = -2.16129
result =
-67/31 275/62 28/31 83/62
======================================================================================================
>>
❤️❤️3. 平方根法
💙函数文件.m
function x=cholesky(A,b)
%喬累斯基分解
[n,n]=size(A);
L=zeros(n,n);%實際上不用為 L 申請空間，使用 A 即可
L(1,1)=sqrt(A(1,1));
for k=2:n
L(k,1)=A(k,1)/L(1,1);
end
for k=2:n-1
L(k,k)=sqrt(A(k,k)-sum(L(k,1:k-1).^2));
for i=k+1:n
L(i,k)=(A(i,k)-sum(L(i,1:k-1).*L(k,1:k-1)))/L(k,k);
end
end
L(n,n)=sqrt(A(n,n)-sum(L(n,1:n-1).^2));
%解下三角方程組 Ly=b
y=zeros(n,1);
for k=1:n
j=1:k-1;
y(k)=(b(k)-L(k,j)*y(j))/L(k,k);
end
%解上三角方程組 L'x=y
x=zeros(n,1);
U=L';
for k=n:-1:1
j=k+1:n;
x(k)=(y(k)-U(k,j)*x(j))/U(k,k);
end
L
U
end
💚测试脚本.m
clear;
clc;
format short;
fprintf('示例：\n');
a = [4 -1 1;-1 4.25 2.75;1 2.75 3.5]
b = [4 6 7.25]'
result = Cholesky(a,b)
fprintf('======================================================================================================\n\n');
💜测试结果
示例：
a =
4.0000 -1.0000 1.0000
-1.0000 4.2500 2.7500
1.0000 2.7500 3.5000
b =
4.0000
6.0000
7.2500
L =
2.0000 0 0
-0.5000 2.0000 0
0.5000 1.5000 1.0000
U =
2.0000 -0.5000 0.5000
0 2.0000 1.5000
0 0 1.0000
result =
1
1
1
======================================================================================================
>>
❤️❤️4. Jacobi迭代
💙函数文件.m
function x=Jacobi(A,b,maxerr,N)
[n,n]=size(A);
x0=zeros(n,1);
x=zeros(n,1); % 给x赋值
k=0;
while k<N
for i=1:n
x(i)=(b(i)-A(i,[1:i-1,i+1:n])*x0([1:i-1,i+1:n]))/A(i,i);
end
if norm(x-x0)<maxerr
break;
end
x0=x;
k=k+1;
disp(['when k=',num2str(k)])
disp('x=');
disp(x); %输出计算的中间结果
end
if k==N
disp('迭代次数已到达上限!');
end
disp(['迭代次数 k=',num2str(k)])
end
💚测试脚本.m
clear;
clc;
format short
fprintf('示例：\n');
a = [2 2 3;4 7 7;-2 4 5]
b = [3 1 -7]'
result = LU(a,b)
fprintf('======================================================================================================\n\n');
💜测试结果
示例：
a =
2 2 3
4 7 7
-2 4 5
b =
3
1
-7
L =
1 0 0
2 1 0
-1 2 1
U =
2 2 3
0 3 1
0 0 6
result =
2
-2
1
======================================================================================================
❤️❤️5. Gauss-Seidel迭代
💙函数文件.m
function x=Gauss_Seidel(A,b,maxerr,N)
[n,n]=size(A);
x0=zeros(n,1);
x=zeros(n,1); % 给x赋值
k=0;
while k<N
for i=1:n
if i==1
x(1)=(b(1)-(A(1,2:n)*x0(2:n,1)))/A(1,1);
elseif i==n
x(n)=(b(1)-(A(n,1:n-1)*x(1:n-1,1)))/A(n,n);
else
x(i)=(b(i)-(A(i,1:i-1)*x(1:i-1,1)+A(i,i+1:n)*x0(i+1:n,1)))/A(i,i);
end
end
if norm(x-x0)<maxerr
break;
end
x0=x;
k=k+1;
disp(['when k=',num2str(k)])
disp('x=');
disp(x); %输出计算的中间结果
end
if k==N
disp('迭代次数已到达上限!');
end
disp(['迭代次数 k=',num2str(k)])
end
💚测试脚本.m
format long
fprintf('示例：\n');
a = [4 -1 1;4 -8 1;-2 1 5]
b = [7 -21 15]'
maxerr = 1e-7;
N=500;
result=Gauss_Seidel(a,b,maxerr,N)
fprintf('=========================================================================================\n\n');
💜测试结果
示例：
a =
4 -1 1
4 -8 1
-2 1 5
b =
7
-21
15
when k=1
x=
1.750000000000000
3.500000000000000
1.400000000000000
when k=2
x=
2.275000000000000
3.937500000000000
1.522500000000000
when k=3
x=
2.353750000000000
3.992187500000000
1.543062500000000
when k=4
x=
2.362281250000000
3.999023437500000
1.545107812500000
when k=5
x=
2.363478906250000
3.999877929687500
1.545415976562500
when k=6
x=
2.363615488281250
3.999984741210938
1.545449247070313
when k=7
x=
2.363633873535156
3.999998092651367
1.545453930883789
when k=8
x=
2.363636040441895
3.999999761581421
1.545454463860474
when k=9
x=
2.363636324430237
3.999999970197678
1.545454535732559
迭代次数 k=9
result =
2.363636358616279
3.999999996274710
1.545454544191570
======================================================================================================
>>
❤️❤️6. 正交多项式曲线拟合
💙函数文件.m
function [] = Curve_fitting(x,y,n)
p = polyfit(x,y,n)
y0=polyval(p,x);
plot(x,y,'o',x,y0,'-') % 作出原数据点 和 拟合曲线图像
% 打印
s='';
len=length(p);
if (len==0)
errordlg('矩阵不能为空');
else
for b=1:len
if b~=len&&b~=len-1
if p(b)<0 && p(b)~=-1
s=strcat(s,num2str(p(b)),'x^',num2str(len-b));
elseif p(b)==-1
s=strcat(s,'-','x^',num2str(len-b));
elseif p(b)>0 && p(b)~=1
if b==1
s=strcat(s,num2str(p(b)),'x^',num2str(len-b));
else
s=strcat(s,'+',num2str(p(b)),'x^',num2str(len-b));
end
elseif p(b)==1
if b==1
s=strcat(s,'x^',num2str(len-b));
else
s=strcat(s,'+','x^',num2str(len-b));
end
else
end
elseif b==len
if p(b)<0 && p(b)~=-1
s=strcat(s,num2str(p(b)));
elseif p(b)==-1
s=strcat(s,'-1');
elseif p(b)>0 && p(b)~=1
s=strcat(s,'+',num2str(p(b)));
elseif p(b)==1
s=strcat(s,'+1');
else
end
else
if p(b)<0 && p(b)~=-1
s=strcat(s,num2str(p(b)),'x');
elseif p(b)==-1
s=strcat(s,'-','x');
elseif p(b)>0 && p(b)~=1
if b==1
s=strcat(s,num2str(p(b)),'x');
else
s=strcat(s,'+',num2str(p(b)),'x');
end
elseif p(b)==1
if b==1
s=strcat(s,'x');
else
s=strcat(s,'+','x');
end
end
end
end
end
fprintf("%d次拟合曲线为：%s\n",n,s)
end
💚测试脚本.m
format short
fprintf('示例：\n');
x = [0 0.25 0.5 0.75 1]
y = [1 1.2840 1.6487 2.117 2.7183]
%x = [5.24 6.15 7.20 7.3 8.9 9.96 12 15.87 19.01 20.87 22.3 24.6 26.7 28.5 29.8 31.35 34.10 35.7]
%y = [3.518 3.264 2.965 2.926 2.542 2.319 1.947 1.403 1.083 0.931 0.815 0.681 0.592 0.512 0.461 0.414 0.34 0.298]
Curve_fitting(x,y,3)
fprintf('======================================================================================================\n\n');
💜测试结果
示例：
x =
0 0.2500 0.5000 0.7500 1.0000
y =
1.0000 1.2840 1.6487 2.1170 2.7183
p =
0.2789 0.4253 1.0141 0.9999
3次拟合曲线为：0.27893x^3+0.42526x^2+1.0141x+0.99991
======================================================================================================
>>
❤️❤️7. 龙贝格积分
💙函数文件.m
function I=Romberg(fun,a,b,e)
k=0; % 迭代次数
n=1; % 区间划分个数
h=b-a; %上下界间距
T(1,1)=h/2*(fun(a)+fun(b));
d=b-a; %误差间距
while e<=d
k=k+1;
h=h/2;
sum=0;
% 计算第一列T
for i=1:n
sum=sum+fun(a+(2*i-1)*h);
end
T(k+1,1)=T(k)/2+h*sum;
% 迭代
for j=1:k
T(k+1,j+1)=T(k+1,j)+(T(k+1,j)-T(k,j))/(4^j-1);
end
n=n*2;
d=abs(T(k+1,k+1)-T(k,k));
end
T
I=T(k+1,k+1);
💚测试脚本.m
clear;
clc;
fprintf('示例：\n');
format long
I=Romberg(@(x) 4/(1+x^2),0,1,0.000001)
fprintf('======================================================================================================\n\n');
💜测试结果
示例：
T =
3.000000000000000 0 0 0 0 0
3.100000000000000 3.133333333333333 0 0 0 0
3.131176470588235 3.141568627450980 3.142117647058823 0 0 0
3.138988494491089 3.141592502458707 3.141594094125888 3.141585783761874 0 0
3.140941612041389 3.141592651224822 3.141592661142563 3.141592638396796 3.141592665277717 0
3.141429893174974 3.141592653552836 3.141592653708037 3.141592653590029 3.141592653649610 3.141592653638244
I =
3.141592653638244
======================================================================================================
>>
——————END-2022-05-22——————
