数值分析中几个简单问题的程序分享

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

和女票一起学习MATLAB，分享几个数值分析问题的实现代码。
q1：Lagrange插值多项式和Newton插值多项式做Runge现象，区间为[-5,5]，将原函数图像和插值函数图象放在同一张图中用不同颜色进行区分。

1. function lagrange(a,b)
   
2. %输入变量:
   
3. %   a:插值多项式的最低次数；
   
4. %   b:插值多项式的最高次数。
   
5. %输出变量：无输出变量，只显示拟合图像与程序运行时间。
   
6. %本函数目的：演示高阶Larange插值的Runge现象。
   
7. %时间：与a,b有关，一般不超过0.1s。
   
8.     %tic;
   
9.     for i=a:b
   
10.       x = linspace(-5, 5, i+1);%
    
11.       y = 1./(1+x.^2);
    
12.       %以上两步确定了插值运算需要的数值对
    
13.       p = polyfit(x,y,i);%p为插值多项式的系数
    
14.       xx = -5:0.01:5;
    
15.       yy = polyval(p,xx);%根据系数计算对应xx的因变量值
    
16.       plot(xx,yy);
    
17.       hold on;
    
18.       grid on;
    
19.     end;
    
20.     plot(xx,1./(1+xx.^2),'r');
    
21.     title('Lagrange插值的Runge现象（红色为实际图像）');
    
22.     %toc;
    
23. end

复制代码

1. function newton(a,b)
   
2. %输入变量:
   
3. %   a:插值多项式的最低次数；
   
4. %   b:插值多项式的最高次数。
   
5. %输出变量：无输出变量，只显示拟合图像与程序运行时间。
   
6. %本函数目的：演示高阶Newton插值的Runge现象。
   
7. %时间：与a,b有关，一般不超过0.3s。
   
8. %备注：需要调用同目录下的calculator_of_newton函数！
   
9.     tic;
   
10.     for i = a:b
    
11.         x = linspace(-5, 5, i+1);
    
12.         y = 1./(1. + x.^2);
    
13.         %以上两步确定了插值运算需要的数值对
    
14.         
    
15.         p = zeros(i+1);
    
16.         p(:,1) = y';%p为差商表与函数值的组合矩阵，第一列为函数值
    
17. 
    
18.         for j1 = 2:i+1
    
19.             for j2 = j1:i+1
    
20.             p(j2, j1) = (p(j2, j1-1) - p(j2-1, j1-1))/(x(j2) - x(j2-j1+1));
    
21.             end
    
22.         end %for循环计算p矩阵完毕
    
23.       
    
24.         x_c = -5:0.01:5;
    
25.         y_c = calculator_of_newton(x_c, p, x);%调用该函数计算自变量经过插值多项式后的对应值
    
26.         plot(x_c, y_c);
    
27.         hold on;
    
28.         grid on;
    
29.     end  
    
30.     plot(x_c, 1./(1. + x_c.^2), 'r');
    
31.     title('Newton插值的Runge现象（红色为实际图像）');
    
32.     toc;
    
33. end

复制代码

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏1 回帖

关键词：简单问题 数值分析 calculator Lagrange function

q2:程序实现复化梯形公式和复化Simpson公式。

1. function y = tra_forml(a, b, n)
   
2. %输入变量;
   
3. %   a:积分下限
   
4. %   b:积分上限
   
5. %   n:划分的区间个数
   
6. %输出变量：
   
7. %   y:数值积分的值
   
8. %本函数的目的：利用复化梯形公式计算给定函数的数值积分
   
9. %时间：不超过0.01s
   
10. %备注：被积函数在本函数同目录下的funcname2.m文件中定义，默认为y = sin(x)/x
    
11.     %tic;
    
12.     x_k = linspace(a, b, n+1);
    
13.     y_k = [];
    
14.     for i = 1:n+1
    
15.         if isnan(funcname2(x_k(i)))
    
16.             syms x;
    
17.             y_k(i) = limit(funcname2(x), x, x_k(i));
    
18.         
    
19.         else
    
20.             y_k(i) = funcname2(x_k(i));
    
21.         end
    
22.     end
    
23.     h = (b - a)/n;
    
24.     y = h./2 * (y_k(1) + 2*sum(y_k(2:n)) + y_k(n+1));
    
25.     %toc;
    
26. end

复制代码

1. function y = spson_forml(a, b, n)
   
2. %输入变量;
   
3. %   a:积分下限
   
4. %   b:积分上限
   
5. %   n:划分的区间个数
   
6. %输出变量：
   
7. %   y:数值积分的值
   
8. %本函数的目的：利用复化Simpson公式计算给定函数的数值积分
   
9. %时间：不超过0.01s
   
10. %备注：被积函数在本函数同目录下的funcname2.m文件中定义，默认为y = sin(x)/x
    
11.     %tic;
    
12.     x_k = linspace(a, b, 2*n+1);
    
13.     y_k1 = [];
    
14.     for i = 1:2:2*n+1
    
15.         if isnan(funcname2(x_k(i)))
    
16.             syms x;
    
17.             y_k1((i+1)./2) = limit(funcname2(x), x, x_k(i));
    
18.         
    
19.         else
    
20.             y_k1((i+1)./2) = funcname2(x_k(i));
    
21.         end
    
22.     end
    
23.     y_k2 = [];
    
24.     for j = 2:2:2*n
    
25.         if isnan(funcname2(x_k(i)))
    
26.             syms x;
    
27.             y_k2(j./2) = limit(funcname2(x), x, x_k(j));
    
28.            
    
29.         else
    
30.             y_k2(j./2) = funcname2(x_k(j));
    
31.         end
    
32.     h = (b - a)./n;
    
33.     y = h./6 * (y_k1(1) + 2*sum(y_k1(2:n)) + 4*sum(y_k2) + y_k1(n+1));
    
34.     %toc;
    
35.     end
    
36. end

复制代码

1. function y = funcname2(x)
   
2. %定义被积函数
   
3. y = sin(x)./x;
   
4. %y = x.^2;
   
5. %y = cos(x);
   
6. %y = (x+1).^3 - 3*x - 4;
   
7. end

复制代码

q3:四阶龙格库塔方法

1. function Runge_Kutta(x0, y0, h, x_end)
   
2. %输入变量：
   
3. %   x0:x的初值
   
4. %   y0:y的初值
   
5. %   h:步长
   
6. %   x_end:x的终止值
   
7. %输出变量：
   
8. %   每一步的（x,y）的值
   
9. %本函数目的：利用四阶Runge_Kutta方法求微分方程的初值问题
   
10. %备注：微分方程在同目录下function3.m中存储，默认为例5.1的问题
    
11. 
    
12.     for i = x0:h:x_end-h %开始循环
    
13.         K1 = function3(x0, y0);
    
14.         K2 = function3(x0 + h/2, y0 + h/2*K1);
    
15.         K3 = function3(x0 + h/2,y0 + h/2*K2);
    
16.         K4 = function3(x0 + h, y0 + h*K3);
    
17.         y1 = y0 + h/6*(K1 + 2*K2 + 2*K3 + K4); %利用K1,K2,K3,K4计算该步长的yn
    
18.         x0 = x0 + h;
    
19.         sprintf('xn = %f, yn = %f',x0, y1) %输出该步的运算结果
    
20.         y0 = y1;
    
21.     end
    
22. end

复制代码

1. function dy = function3(x0, y0)
   
2.     dy = y0 - 2*x0/y0;
   
3. end

复制代码

q4:方程求根（二分法，牛顿法）

1. function x_med = dichotomy(a, b, tol)
   
2. %输入变量;
   
3. %   a:区间下限
   
4. %   b:区间上限
   
5. %   tol:可选参数，用于控制误差，默认为1e-6
   
6. %输出变量：
   
7. %   x_med:求得的方程的根
   
8. %另外输出：
   
9. %   迭代次数
   
10. %本函数的目的：利用二分法解方程
    
11. %时间：不超过0.001s
    
12. %备注：被积函数在本函数同目录下的funcname4.m文件中定义
    
13.     %tic;
    
14.     if nargin == 2
    
15.         %根据传入参数个数，判断是否有tol传入
    
16.         %nargin = 2，说明只传入a, b的值，此时设置tol的默认值
    
17.         tol = 1e-6;
    
18.     end
    
19.     x1 = a;
    
20.     x2 = b;
    
21.     y1 = funcname4(x1);
    
22.     y2 = funcname4(x2);
    
23.    
    
24.     if y1 * y2 > 0 %若区间断点内的所有点（包括边界）的函数值同号，则根据单调性区间内无解
    
25.         fprintf('本区间内无解！\n');
    
26.         return;
    
27.     else
    
28.         x_med = (x1 + x2)/2;
    
29.         y_med = funcname4(x_med);
    
30.         N = 0;%计算迭代次数，每次执行while则+1
    
31.         while abs(y_med) > tol %若区间中间值的函数值与0的距离小于给定误差，则终止
    
32.             x_med = (x1 + x2)/2;
    
33.             y_med = funcname4(x_med);
    
34.         
    
35.             if y2 * y_med > 0
    
36.                 x2 = x_med;
    
37.                 y2 = y_med;
    
38.             else
    
39.                 x1 = x_med;
    
40.                 y1 = y_med;
    
41.             end
    
42.             N = N + 1;
    
43.         end
    
44.     end
    
45.     fprintf('迭代次数为');
    
46.     disp(N)
    
47.     %toc;
    
48. end

复制代码

1. function x_f = newton_solver(x0, show, tol)
   
2. %输入变量;
   
3. %   x0:方程的近似初始解
   
4. %   show:是否显示图像逼近过程，默认为1，表示显示，可传入其他参数关闭图像展示
   
5. %   tol:可选参数，用于控制误差，默认为1e-6
   
6. %输出变量：
   
7. %   x_f:求得的方程的根
   
8. %另外输出：   
   
9. %   展示逼近过程的图像（if show == 1）
   
10. %   迭代次数
    
11. %本函数的目的：利用newton法解方程
    
12. %时间：不超过0.001s
    
13. %备注：被积函数在本函数同目录下的funcname4.m文件中定义,
    
14. %      为演示效果明显,采用funcname.m中的y=(x+1).^3-3*x-4函数,初始值在[2,3)内选择
    
15.     if nargin == 1
    
16.         %根据传入参数个数，判断是否采用默认值
    
17.         %nargin == 1，说明只传入x0的值，此时设置show, tol的默认值
    
18.         show = 1;
    
19.         tol = 1e-6;
    
20.     end
    
21.    
    
22.     if nargin == 2
    
23.         %nargin == 2，说明只传入x0,show的值，此时设置 tol的默认值
    
24.         tol = 1e-6;
    
25.     end
    
26. 
    
27.     f0 = funcname(x0);
    
28.     syms s;
    
29.     df = diff(funcname(s));%求出待求解函数的导函数
    
30.    
    
31.     if abs(f0) < tol
    
32.         %若初始值的函数值小于给定误差限，则直接返回初始值作为根
    
33.         x_f = x0;
    
34.         return;
    
35.     else
    
36.         N = 0;%计算迭代次数，每次执行while则+1
    
37.         while abs(f0) > tol %当求得函数值与0的距离小于给定误差时，终止
    
38.             df0 = subs(df, s, x0);
    
39.             b0 = f0 - df0*x0;
    
40.             if df0 == 0
    
41.                 fprintf('在该初始值处不收敛，请更换初始值\n');
    
42.                 return;
    
43.             end
    
44.             x1 = x0 - f0/df0;
    
45.             f1 = funcname(x1);
    
46.             if show == 1 %若打开图像显示过程
    
47.                 %画出原始方程图像
    
48.                 axis([0.5 3 0 50]);%设置画布
    
49.                 x = 0:0.01:3;
    
50.                 y = funcname(x);
    
51.                 plot(x,y);%画图
    
52.                 xlabel('X');ylabel('Y');ylim([0 50]);
    
53.                 annotation('arrow',[0.132 0.132],[0.8 1]);%添加坐标轴箭头
    
54.                 annotation('arrow',[0.8 1],[0.108 0.108]);%添加坐标轴箭头
    
55.                 hold on;
    
56.    
    
57.                 plot(x0, f0, '.','MarkerEdgeColor','r','MarkerSize',16);%画出初始点的位置
    
58.                 pause(0.5);%暂停0.5s
    
59.                 s1line = refline(double(df0), double(b0));%添加切线
    
60.                 set(s1line,'Color','r');ylim([0 50]);xlim([0.5 3])%设置切线显示范围
    
61.    
    
62.                 plot(x1, 0, '.', 'MarkerEdgeColor','r','MarkerSize',16); %画出切线与x轴交点
    
63.                 pause(0.5);
    
64.                 line([x1 x1], [0 50]);%添加竖直线
    
65.                 plot(x1,f1, '.', 'MarkerEdgeColor','r','MarkerSize',16);%画出竖直线与原函数的交点
    
66.                 pause(0.5);
    
67.                 hold off;
    
68.                 pause(2);
    
69.             end
    
70.             %以下为迭代过程
    
71.             x0 = x1;
    
72.             f0 = f1;
    
73.             df0 = subs(df, s, x0);
    
74.             x1 = x0 - f0/df0;
    
75.             f1 = funcname(x1);
    
76.             N = N + 1;
    
77.         end
    
78.         x_f = double(x0);
    
79.     end
    
80.     fprintf('迭代次数为');
    
81.     disp(N)
    
82. end

复制代码

1. function y = funcname4(x)
   
2. %定义待求解方程函数
   
3. y = (x+1).^3 - 3*x - 4;
   
4. end

复制代码

q5:解线性方程
1.LU分解

1. function [L,U] = LU_deco(A)
   
2. %输入变量：
   
3. %   A:非奇异系数矩阵
   
4. %输出变量：
   
5. %   L:下三角矩阵
   
6. %   U:上三角矩阵
   
7. %本函数的目的：对线性方程组的非奇异系数矩阵进行LU分解
   
8. %
   
9.     [n n] = size(A);%获得矩阵A的维度
   
10.     U = A;%初始化矩阵U
    
11.     L = eye(n);%初始化矩阵L
    
12.     for k = 1:n-1 %k控制消元次数
    
13.         for raw = (k+1):n %raw为当前消元次数下需要更新的行
    
14.             L(raw,k) =  U(raw,k)/U(k,k); %矩阵重新赋值
    
15.             U(raw,:) = U(raw,:) - U(raw,k)/U(k,k).*U(k,:); %矩阵重新赋值
    
16.         end
    
17.     end
    
18. end

复制代码

2.1 Jacobi迭代法求解

1. function x = Jacobi(A,b,k)
   
2. %输入变量：
   
3. %   A:线性方程组的系数矩阵
   
4. %   b:常向量
   
5. %   k:设置的迭代次数
   
6. %输出变量：
   
7. %   迭代计算后的解
   
8. %本函数目的：利用Jacobi迭代法解线性方程组
   
9.     [n n] = size(A); %获取未知数个数
   
10.     x = rand(1,n); %随机化初始解
    
11.     i = 1; %循环因子
    
12.     x_ = zeros(1,n); %用于存储迭代一次后的新的解
    
13.     while i <= k %迭代开始
    
14.         for j =1:n
    
15.             x_(j) = (b(j) - sum(A(j,:).*x) + A(j,j)*x(j))/A(j,j);         
    
16.         end
    
17.         x = x_; %用迭代后的解重新赋值于x
    
18.         i = i + 1;
    
19.         if mod(i,3) == 0 %控制输出
    
20.             sprintf('第 %d 次迭代后的解为：\n',i)
    
21.             disp(x_);
    
22.         end      
    
23.     end
    
24. end

复制代码

2.2 Gauss-Seidel迭代法求解

1. function x = Gauss_Seidel(A,b,k)
   
2. %输入变量：
   
3. %   A:线性方程组的系数矩阵
   
4. %   b:常向量
   
5. %   k:设置的迭代次数
   
6. %输出变量：
   
7. %   迭代计算后的解
   
8. %本函数目的：利用Gauss-Seidel迭代法解线性方程组
   
9.     [n n] = size(A); %获取未知数个数
   
10.     x = rand(1,n); %随机化初始解
    
11.     i = 1; %循环因子
    
12.     x_ = zeros(1,n); %用于存储迭代一次后的新的解
    
13.     while i <= k %迭代开始
    
14.         for j =1:n
    
15.             x_(j) = (b(j) - sum(A(j,:).*x) + A(j,j)*x(j))/A(j,j);
    
16.             x = x_; %解向量的下一个值均用最新迭代后的值求解
    
17.         end
    
18.         
    
19.         i = i + 1;
    
20.         if mod(i,3) == 0 %控制输出
    
21.             sprintf('第 %d 次迭代后的解为：\n',i)
    
22.             disp(x_);
    
23.         end      
    
24.     end
    
25. end

复制代码

q6:利用幂法求解矩阵的最大特征值

1. function [Eigenvalue Eigenvector] = Eigenvalues(A, k)
   
2. %输入变量：
   
3. %   A:待求解的矩阵
   
4. %   k:设置的迭代次数
   
5. %输出变量：
   
6. %   Eigenvalue：最大特征值
   
7. %   Eigenvector：对应Eigenvalue的特征向量
   
8. %本函数目的：利用幂法求解矩阵的最大特征值
   
9.     [n n] = size(A); %获取矩阵维度
   
10.     v = rand(n,1); %初始化向量，v0
    
11.     i = 1; %迭代控制变量
    
12.     v = A * v;  %v1
    
13.     while i <= k %迭代开始
    
14.         u = v/max(abs(v));
    
15.         v = A*u;
    
16.         i = i+1;
    
17.     end
    
18.     Eigenvalue = max(v); %返回值
    
19.     Eigenvector = u;
    
20. end

复制代码

总结一下：这是本人第一次系统学习MATLAB，深切体会到了各种软件在各自所擅长的领域为计算带来的方便！另外上述代码只是在一些相对较苛刻的条件满足下才能正确运行，值得改进的地方很多，但作为新手参考仍有价值。另外也可以解决同一问题的不同实现代码，甚至MATLAB中的源码，这样可以加速理解。
最后，理论知识很重要！理论知识很重要！理论知识，真的很重要！

您好！谢谢您的程序，很受益！不知是否可将calculator_of_newton.m发到yyf_brave@163.com，谢谢