Matlab学习笔记（三）

主要内容（大概）

一、复习积分、级数
二、泰勒级数和傅里叶级数
三、Matlab求解 —— 微分和代数

目录

积分复习

级数复习

泰勒级数

傅里叶级数

Matlab求解微分方程

积分复习 ↶

1.quad做无穷（广义）有点问题，int比较精确，用int最好，下面的例子用来对比：
例，计算下列广义积分
$${\int }_{-\infty }^{\infty }{e}^{cosx-x^2}dx$$

clear
format long
syms x %为int做准备
ff = @(x)(exp(cos(x)-x.^2)); %定义匿名函数，可以用inline
y1 = quad(ff, -inf, inf)
y2 = quad(ff, -100, 100)
y3 = quad(ff, -1e10, 1e10)
y4 = quadl(ff, -1e20, 1e20)
y5 = vpa(int(exp(cos(x)-x.^2), -inf, inf), 15)

y1 =  NaN
y2 = 3.989813604029027
y3 = 3.989812541283583
y4 = 25.814543367399441 %精度的范围导致的结果变化
y5 = 3.989812272586 %用vpa控制精度产生的结果

2.quad/quadl（精度问题）优化程序
优化程序：

clear;
ff = @(x)(exp(cos(x)-x.^2));
n = 10; m = 2; c = 1e-6; a = inf; % n为上下限，m控制步长，c指定精度，a为辅助变量
b = quadl(ff, -n, n);
while abs(a-b) > c;
a = b; n = n*m;
b = quadl(ff, -n, n);
end

级数复习 ↶

应用一：Euler常数的发现
$$\gamma =\underset{n\to \infty }{\lim }(\sum _{i=1}^n\frac{1}{i}-ln(n))$$
1.观察调和级数部分和数列的规律

clear;
syms m
for n=1:100;
a(n) = symsum(1/m, m, 1, n);
end
plot([1:n], a, 'o')

2.画图比较与对数函数的差异

hold on
for n=1:100;
b(n) = log(n);
end
plot([1:n], b, '*')

泰勒级数 ↶

应用二：Taylor幂级数展开
例：计算cos(2x)泰勒x=4处展开式的前6项
第一步：编辑Ty.m函数文件

function y = Ty(fun, a, n)
syms x
c = subs(fun, x, a); %用a替换x
for k=1:n
funk = subs(diff(fun, k), x, a);
c = [c, funk]; %有加法的感觉，把funk放入c中， c = 0; funk = 1; c += funk
end
d = 1;
for k=1:n
d = [d, (x-a)^k/factorial(k)]; %同理
end
md = c.*d; %两者相乘
y = sum(md, 2); %对md的每行分别求和

第二步：命令窗口输入

syms x
fun = cos(2*x)
Ty(fun, 4, 5)

结果：

fun = cos(2*x)
ans = cos(8) + (4*sin(8)*(x - 4)^3)/3 - (4*sin(8)*(x - 4)^5)/15 - 2*sin(8)*(x - 4) - 2*cos(8)*(x - 4)^2 + (2*cos(8)*(x - 4)^4)/3

Talylor幂级数展开

syms x
taylor(fun, x, 'ExpansionPoint', x0, 'order', n)
%或者
taylor(fun, x, x0, 'order', n)
%幂级数
taylor(fun, x, n)

例：$$f(x)=\frac{sinx}{x^2+4x+3}$$，求其Taylor幂级数展开前9项，并求关于x=2和x=a进行Taylor幂级数展开

syms x;
f = sin(x)/(x^2+4*x+3);
y1 = taylor(f, x, 9) 
%结果
y1 = sin(9)/120 + (x - 9)*(cos(9)/120 - (11*sin(9))/7200) - (x - 9)^2*((11*cos(9))/7200 + (1709*sin(9))/432000) - (x - 9)^3*((509*cos(9))/432000 - (19129*sin(9))/25920000) + ((58051*cos(9))/1555200000 - (4763231*sin(9))/93312000000)*(x - 9)^5 + ((5929*cos(9))/25920000 + (380851*sin(9))/1555200000)*(x - 9)^4

x = 2时：

syms x;
f = sin(x)/(x^2+4*x+3);
y1 = taylor(f, x, 2, 'order', 9)
%结果
y1 = sin(2)/15 - ((131623*cos(2))/35880468750 + (875225059*sin(2))/34445250000000)*(x - 2)^8 + (x - 2)*(cos(2)/15 - (8*sin(2))/225) - ((623*cos(2))/11390625 + (585671*sin(2))/2733750000)*(x - 2)^6 + ((262453*cos(2))/19136250000 + (397361*sin(2))/5125781250)*(x - 2)^7 - (x - 2)^2*((8*cos(2))/225 + (127*sin(2))/6750) + (x - 2)^3*((23*cos(2))/6750 + (628*sin(2))/50625) + (x - 2)^4*((28*cos(2))/50625 - (15697*sin(2))/6075000) + (x - 2)^5*((203*cos(2))/6075000 + (6277*sin(2))/11390625)

x = a时；

syms x a;
f = sin(x)/(x^2+4*x+3);
y1 = taylor(f, x, a, 'order', 9)
%结果
y1 = sin(a)/(a^2 + 4*a + 3) - (a - x)^2*(sin(a)/(2*(a^2 + 4*a + 3)) + (cos(a)*(2*a + 4))/(a^2 + 4*a + 3)^2 - (sin(a)*((2*a + 4)^2/(a^2 + 4*a + 3)^2 - 1/(a^2 + 4*a + 3)))/(a^2 + 4*a + 3)) - (cos(a)/(a^2 + 4*a + 3) - (sin(a)*(2*a + 4))/(a^2 + 4*a + 3)^2)*(a - x) - (a - x)^6*(sin(a)/(720*(a^2 + 4*a + 3)) - (sin(a)*(((2*a + 4)^2/(a^2 + 4*a + 3)^2 - 1/(a^2 + 4*a + 3))/(a^2 + 4*a + 3) + ((2*a + 4)*((2*a + 4)/(a^2 + 4*a + 3)^2 - (((2*a + 4)^2/(a^2 + 4*a + 3)^2 - 1/(a^2 + 4*a + 3))*(2*a + 4))/(a^2 + 4*a + 3)))/(a^2 + 4*a + 3)))/(2*(a^2 + 4*a + 3)) + (cos(a)*((2*a + 4)/(a^2 + 4*a + 3)^2 - (((2*a + 4)^2/(a^2 + 4*a + 3)^2 - 1/(a^2 + 4*a + 3))*(2*a + 4))/(a^2 + 4*a + 3)))/(6*(a^2 + 4*a + 3)) + (cos(a)*(2*a + 4))/(120*(a^2 + 4*a + 3)^2) - (sin(a)*((2*a + 4)^2/(a^2 + 4*a + 3)^2 - 1/(a^2 + 4*a + 3)))/(24*(a^2 + 4*a + 3)) - (sin(a)*((((2*a + 4)^2/(a^2 + 4*a + 3)^2 - 1/(a^2 + 4*a + 3))/(a^2 + 4*a + 3) + ((2*a + 4)*((2*a + 4)/(a^2 + 4*a + 3)^2 - (((2*a + 4)^2/(a^2 + 4*a + 3)^2 - 1/(a^2 + 4*a + 3))*(2*a + 4))/(a^2 + 4*a + 3)))/(a^2 + 4*a + 3))/(a^2 + 4*a + 3) + ((2*a + 4)*(((2*a + 4)/(a^2 + 4*a + 3)^2 - (((2*a + 4)^2/(a^2 + 4*a + 3)^2 - 1/(a^2 + 4*a + 3))*(2*a + 4))/(a^2 + 4*a + 3))/(a^2 + 4*a + 3) - ((2*a + 4)*(((2*a + 4)^2/(a^2 + 4*a + 3)^2 - 1/(a^2 + 4*a + 3))/(a^2 + 4*a + 3) + ((2*a + 4)*((2*a + 4)/(a^2 + 4*a + 3)^2 - (((2*a + 4)^2/(a^2 + 4*a + 3)^2 - 1/(a^2 + 4*a + 3))*(2*a + 4))/(a^2 + 4*a + 3)))/(a^2 + 4*a + 3)))/(a^2 + 4*a + 3)))/(a^2 + 4*a + 3)))/(a^2 + 4*a + 3) + (cos(a)*(((2*a + 4)/(a^2 + 4*a + 3)^2 - (((2*a + 4)^2/(a^2 + 4*a + 3)^