当前位置:   article > 正文

数值分析——Gauss-Legendre 求积分(Matlab实现)_matlab中有gauss-legendre函数吗

matlab中有gauss-legendre函数吗

2020 4.4

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

∫ 0 1 sin ⁡ x x d x \int_{0}^{1} \frac{\sin x}{x} d x 01xsinxdx

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

n x k A k n x k A k 0 0 2 3 ± 0.8611363116 0.3478548451 1 ± 0.5773502692 1 ± 0.3399810436 0.6521451549 2 ± 0.7745966692 0.5555555555 4 ± 0.9061798459 0.2369268851 0 0.8888888888 ± 0.5384693101 0.4786286705 0 0.568888889 nxkAknxkAk0023±0.86113631160.34785484511±0.57735026921±0.33998104360.65214515492±0.77459666920.55555555554±0.90617984590.236926885100.8888888888±0.53846931010.478628670500.568888889 n012xk0±0.5773502692±0.77459666920Ak210.55555555550.8888888888n34xk±0.8611363116±0.3399810436±0.9061798459±0.53846931010Ak0.34785484510.65214515490.23692688510.47862867050.568888889

区间变换与函数变换:

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

函数变换:
∫ a b f ( x ) d x = b − a 2 ∫ − 1 1 f [ 1 2 ( a + b ) + 1 2 ( b − a ) t ] d 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 abf(x)dx=2ba11f[21(a+b)+21(ba)t]dt

具体实现

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


  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29

踩坑日志

错误:

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

出错 sym/subsref (line 870)
            R_tilde = builtin('subsref',L_tilde,Idx);
  • 1
  • 2
  • 3
  • 4
  • 5

原因:

I=2*f1(0);
  • 1

应该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;
  • 1
  • 2
  • 3
  • 4
  • 5

原因:

syms x
  • 1

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

改进:

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

这样可以不用声明变量

学到的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
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24

matlab函数的编写

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

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

函数语句的语法是:

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

例:

建立函数文件,命名为 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
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16

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

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

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

help mymax
  • 1

MATLAB执行上述语句,返回以下结果:

This function calculates the maximum of the
 five numbers given as input
  • 1
  • 2

可以调用该函数:

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

MATLAB执行上述语句,返回以下结果:

ans =
    89
  • 1
  • 2

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

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

I =

   0.946041135536224

  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11

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

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

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

元胞数组

声明:本文内容由网友自发贡献,不代表【wpsshop博客】立场,版权归原作者所有,本站不承担相应法律责任。如您发现有侵权的内容,请联系我们。转载请注明出处:https://www.wpsshop.cn/w/知新_RL/article/detail/811272
推荐阅读
相关标签
  

闽ICP备14008679号