数值分析——Gauss-Legendre 求积分(Matlab实现)

2020 4.4

题目:编写Gauss求积法计算积分的程序（Gauss点数取1,2,3,4,5即可）并用于计算积分

$${\int }_0^1\frac{\sin x}{x}dx$$

部分高斯-勒让德求积公式的结点和求积系数:

$$\begin{array}{cccccc}n& x_k& A_k& n& x_k& A_k\\ 0& 0& 2& 3& \pm 0.8611363116& 0.3478548451\\ 1& \pm 0.5773502692& 1& & \pm 0.3399810436& 0.6521451549\\ 2& \pm 0.7745966692& 0.5555555555& 4& \pm 0.9061798459& 0.2369268851\\ & 0& 0.8888888888& & \pm 0.5384693101& 0.4786286705\\ & & & & 0& 0.568888889\end{array}$$

区间变换与函数变换:

区间[a,b]->[-1,1]:
$$x=\frac{1}{2}(a+b)+\frac{1}{2}(b-a)t$$

函数变换:
$${\int }_a^{b}f(x)\mathrm{d}x=\frac{b-a}{2}{\int }_{-1}^{1}f\left[\frac{1}{2}(a+b)+\frac{1}{2}(b-a)t\right]\mathrm{d}t$$

具体实现

function [f1,I] = GaussInt(f,a,b,n)
%该函数利用Gauss-Lagendre求积公式计算积分
%依次输入函数,左右区间端点,和结点个数进行计算高斯积分
syms t
f1=subs(f,symvar(f),(b+a)/2+(b-a)*t/2)*(b-a)/2;
switch n %分别代表不同结点个数时的Gauss积分
    case 1
    I=eval(2*subs(f1,symvar(f1),0));
    case 2
    I=eval(subs(f1,symvar(f1),0.5773503))+...
      eval(subs(f1,symvar(f1),-0.5773503));
    case 3
    I=0.5555555555*eval(subs(f1,symvar(f1),0.7745966692))+...
      0.5555555555*eval(subs(f1,symvar(f1),-0.7745966692))+...
      0.8888888888*eval(subs(f1,symvar(f1),0));
    case 4
    I=0.3478548451*eval(subs(f1,symvar(f1),0.8611363116))+...
      0.3478548451*eval(subs(f1,symvar(f1),-0.8611363116))+...
      0.6521451549*eval(subs(f1,symvar(f1),0.3399810436))+...
      0.6521451549*eval(subs(f1,symvar(f1),-0.3399810436));
    case 5
     I=0.2369268851*eval(subs(f1,symvar(f1),0.9061798459))+...
      0.2369268851*eval(subs(f1,symvar(f1),-0.9061798459))+...
      0.4786286705*eval(subs(f1,symvar(f1),0.5384693101))+...
      0.4786286705*eval(subs(f1,symvar(f1),-0.5384693101))+...
      0.568888889*eval(subs(f1,symvar(f1),0));
end

踩坑日志

错误:

>> [F,I]=GaussInt(sin(x)/x, 0, 1,1)
数组索引必须为正整数或逻辑值。

出错 sym/subsref (line 870)
            R_tilde = builtin('subsref',L_tilde,Idx);

原因:

I=2*f1(0);

应该f1是一个符号表达式而不是一个函数

错误:

>> [F,I]=GaussInt(sin(x)/x, 0, 1,1)
未定义函数或变量 'x'。

出错 GaussInt (line 5)
f1=subs(f,x,(b-a)/2+(b-a)*t/2)*(b-a)/2;

原因:

syms x

所编写的函数中没有定义变量x

改进:

f1=subs(f,symvar(f),(b-a)/2+(b-a)*t/2)*(b-a)/2;

这样可以不用声明变量

学到的matlab语法

matlab代码注释

%

快速:

Ctrl +t

Ctrl+ r

switch多分支结构

switch n %分别代表不同结点个数时的Gauss积分
    case 1
    I=eval(2*subs(f1,symvar(f1),0));
    case 2
    I=eval(subs(f1,symvar(f1),0.5773503))+...
      eval(subs(f1,symvar(f1),-0.5773503));
    case 3
    I=0.5555555555*eval(subs(f1,symvar(f1),0.7745966692))+...
      0.5555555555*eval(subs(f1,symvar(f1),-0.7745966692))+...
      0.8888888888*eval(subs(f1,symvar(f1),0));
    case 4
    I=0.3478548451*eval(subs(f1,symvar(f1),0.8611363116))+...
      0.3478548451*eval(subs(f1,symvar(f1),-0.8611363116))+...
      0.6521451549*eval(subs(f1,symvar(f1),0.3399810436))+...
      0.6521451549*eval(subs(f1,symvar(f1),-0.3399810436));
    case 5
     I=0.2369268851*eval(subs(f1,symvar(f1),0.9061798459))+...
      0.2369268851*eval(subs(f1,symvar(f1),-0.9061798459))+...
      0.4786286705*eval(subs(f1,symvar(f1),0.5384693101))+...
      0.4786286705*eval(subs(f1,symvar(f1),-0.5384693101))+...
      0.568888889*eval(subs(f1,symvar(f1),0));
    otherwise
   	  disp("结点数为1-5")
end

matlab函数的编写

在MATLAB中，函数定义在单独的文件。文件函数的文件名应该是相同的。

函数可以接受多个输入参数和可能返回多个输出参数。

函数语句的语法是：

function [out1,out2, ..., outN] = myfun(in1,in2,in3, ..., inN)

例:

建立函数文件，命名为 mymax.m并输入下面的代码：

function max = mymax(n1, n2, n3, n4, n5)
%This function calculates the maximum of the
% five numbers given as input
max =  n1;
if(n2 > max)
    max = n2;
end
if(n3 > max)
   max = n3;
end
if(n4 > max)
    max = n4;
end
if(n5 > max)
    max = n5;
end

每个函数的第一行要以 function 关键字开始。它给出了函数的名称和参数的顺序。

在我们的例子中，mymax 函数有5个输入参数和一个输出参数。

注释行语句的功能后提供的帮助文本。这些线条打印，当输入：

help mymax

MATLAB执行上述语句，返回以下结果：

This function calculates the maximum of the
 five numbers given as input

可以调用该函数：

mymax(34, 78, 89, 23, 11)

MATLAB执行上述语句，返回以下结果：

ans =
    89

注:如果要输出多个返回值:

>> [F,I]=GaussInt(sin(x)/x, 0, 1,2)
 
F =
 
sin(t/2 + 1/2)/(2*(t/2 + 1/2))
 

I =

   0.946041135536224

一种matlab的数据类型能够存储变长的一维向量:

刚开始想把结点和系数存储起来用循环进行对应的遍历

A=[1,2,3,4];
B=[5];
C=[6,7];
p={A,B,C}
以下运行结果：
>> p{1}(1,2)
ans =
2

元胞数组