教了个数值分析
本学期我教了个数值分析和matlab课程，当然matlab很久没有碰过了，这次是临时重学了一遍。至于数值分析，我是一点也不怕的，因为矩阵论部分内容刚刚教过。当然，如果不深入讲解它的理论，知识告诉你怎么去用它，我想这些就绰绰有余了。
重点：
1.范数
2.线性方程组求解的迭代法
3.LU分解，乔列斯基分解，三对角方程组
4.Lagrang插值和Newton插值，样条插值。
5.数值积分的计算 以前印象深刻的是广义皮亚诺定理， 当然现在发现用Euler-Maclarin公式也是很有意思的事情，这对于做积分估计是很有用的，同时广义皮亚诺定理的证明实际就 是选择Hermit插值，特别的对于Simpson公式的证明，硬用大一微积分的做法可能说不清楚，但是引入带导数的插值后，一切迎刃而解。
6.非线性方程求解的迭代法，主要是一元清醒。事实上，多元情形也是对的，包括所谓的Jocobi，Gauss-Seidel迭代也可推广至非线性方程组情形。在多元情形需要利用隐函数定理。
7.常微分方程的数值解法
大概就这些吧
下面是我自己编写的一些程序，当然我是尽量结合C语言的风格来编写的，当然matlab在赋值和变量类型方面有很多比C语言更直观的处理，接近数学运算。
顺序高斯消元法的matlab程序: m文件
%Gauss.m
function [x,A,b]=Gauss(A,b)
n=length(b);
% 顺序消元法
for k=1:(n-1)
for i=(k+1):1:n
b(i)=b(i)-b(k)*(A(i,k)/A(k,k));
for j=n:-1:k
A(i,j)=A(i,j)-A(k,j)*(A(i,k)/A(k,k));
end
end
end
%回代
x=zeros(n,1);
x(n)=b(n)/A(n,n);
for k=n-1:-1:1
s=0;
for m=(k+1):n
s=s+A(k,m)*x(m);
end
x(k)=(b(k)-s)/A(k,k);
end
以下在命令窗口调用:
>> A=[2 -1 1; 1 1 1; -1 -1 2];
>> b=[2 3 0]';
>> Gauss(A,b)
ans =
1
1
1
Jacobi迭代法的矩阵形式matlab程序，这是在别的书上参考的。
function [x,number]=Jacobidiedai(A,b,x0,er)
%Jacobi迭代法, x0是初值, er是误差, number是迭代次数；
D=diag(diag(A));
DD=inv(D);
U=triu(A,1);
L=tril(A,-1);
B=-DD*(L+U);
f=DD*b;
number=0;
x=B*x0+f;
while norm(x-x0)>er
x0=x;
x=B*x0+f;
number=number+1;
end
在命令窗口调用
>> A=[5,1,1;2,7,1;1,2,9];
>> b=[7,10,12]';
>> x0=[0,0,0]';
>> er=1e-6;
>> [x,number]=Jacobidiedai(A,b,x0,er)
x =
1.00000011158925
1.00000011819469
1.00000009922229
number =
16
求解Lagrange插值公式的matlab程序 m文件
% lagpoly.m
%求解Lagrange插值多项式
function lag=lagploy(x1,y1)
% x1是插值节点, y1是插值节点处的函数值；
n=length(x1);
p=zeros(1,n);
for i=1:n
t=1;
for j=1:n
if j~=i
t=t/(x1(i)-x1(j));
end
end
p(i)=y1(i)*t;
end
syms x
q=0;
for i=1:n
s=1;
for j=1:n
if j~=i
s=s*(x-x1(j));
end
end
q=q+p(i)*s;
end
lag=simplify(q);
在命令窗口调用
>> x=[1 2 3 4];
>> y=[1 8 27 64];
>> lagpoly(x,y)
ans =
x^3
Newton法解方程的m文件
function [y,k]=egnewton(x0,er)
f=inline('x^3+x-1');
syms x
g=diff(f(x));
y(1)=x0;
x1=x0-f(x0)/(subs(g,x,x0));
y(2)=x1;
k=2;
while abs(x1-x0)>er
x1=x0-f(x0)/(subs(g,x,x0));
x2=x1-f(x1)/(subs(g,x,x1));
k=k+1;
y(k)=x2;
x0=x1;
x1=x2;
end
命令窗口调用
>> egnewton(0.5,1e-6)
ans =
0.50000000000000 0.71428571428571 0.68317972350230 0.68232842330458 0.68232780382835
>> [y,k]=egnewton(0.5,1e-6)
y =
0.50000000000000 0.71428571428571 0.68317972350230 0.68232842330458 0.68232780382835
k =
5
将上述m文件改变一点点更加整齐
function [y,k]=haha(x0,er)
syms x
f=x^2-2;
g=diff(f);
y(1)=x0;
x1=x0-subs(f,x,x0)/(subs(g,x,x0));
y(2)=x1;
k=2;
while abs(x1-x0)>er
x1=x0-subs(f,x,x0)/(subs(g,x,x0));
x2=x1-subs(f,x,x1)/(subs(g,x,x1));
k=k+1;
y(k)=x2;
x0=x1;
x1=x2;
end
haha(1,0.0001)
ans =
1.00000000000000 1.50000000000000 1.41666666666667 1.41421568627451 1.41421356237469
二分法求解
function y=bicut(m,n,er)
a=m;
b=n;
k=0;
ff=inline('x^3+x-1');
while abs(b-a)>er
xk=(b+a)/2;
fk=ff(xk);
fa=ff(a);
k=k+1;
if fk==0
y(k)=xk;
break;
elseif fa*fk<0
b=xk;
else
a=xk;
end
y(k)=xk;
end
bicut(-1,2,1e-6)
ans =
Columns 1 through 5
0.50000000000000 1.25000000000000 0.87500000000000 0.68750000000000 0.59375000000000
Columns 6 through 10
0.64062500000000 0.66406250000000 0.67578125000000 0.68164062500000 0.68457031250000
Columns 11 through 15
0.68310546875000 0.68237304687500 0.68200683593750 0.68218994140625 0.68228149414063
Columns 16 through 20
0.68232727050781 0.68235015869141 0.68233871459961 0.68233299255371 0.68233013153076
Columns 21 through 22
0.68232870101929 0.68232798576355
复化求积公式
function [s,er]=complexint(t,a,b,n)
dx=(b-a)/n;
syms x
ff=4/(1+x^2);
xx=a:dx:b;
s=0;
if t==1 % 复化梯形求积公式情形
for i=1:n
f1=subs(ff,x,xx(i)); %小区间左端点的函数值
f2=subs(ff,x,xx(i+1)); %小区间右端点的函数值
s=s+(f1+f2)*dx/2;
end
end
if t==2 % 复化梯形求积公式情形
for i=1:n
f1=subs(ff,x,xx(i)); %小区间左端点的函数值
f2=subs(ff,x,xx(i+1)); %小区间右端点的函数值
f12=subs(ff,x,(xx(i)+xx(i+1))/2); %小区间中点的函数值
s=s+(f1+4*f12+f2)*dx/6;
end
end
er=abs(pi-s); % 计算结果与真实值的误差
>> [s,er]=complexint(1,0,1,5)
s =
3.13492611381099
er =
0.00666653977880
>> [s,er]=complexint(2,0,1,5)
s =
3.14159261393922
er =
3.965057793209326e-008
编写了一个约瑟夫问题的matlab程序，昨天看日记知道2004年时自己用数组编写了一个C程序，但是没有记录下来，现在用matlab编写了一个，当然还是用数组的想法，自然还是模拟一遍。显然这个程序不是最优的。建立m文件函数如下：
function z=Josephus(n,m);
x=ones(1,n);
y=zeros(1,n);
i=1;
j=0;
while isequal(x,y)~=1
if i<n
if (x(i))~=0
j=j+1;
end
i=i+1;
if (rem(j+1,m)==0 & x(i)~=0)
x(i)=0;
disp(i)
j=0;
if i<n
i=i+1;
else
i=1;
end
end
else
if (x(i))~=0
j=j+1;
end
i=1;
if (rem(j+1,m)==0 & x(1)~=0)
x(1)=0;
disp(1)
j=0;
i=2;
end
end
end
在命令窗口调用 Josephus(8,3) 可得
>> Josephus(8,3)
3
6
1
5
2
8
4
7
以上给出了被铲除的序号，最后一个被铲除的自然是最后剩下的了，即7号是所要求的数字。
这倒使我想起一件事情，以前数学建模时有个叫王青的讲解了隔2舍的一个递推公式，划归为2进制，而这是baidu上也有的算法，而一般情形，显然baidu上的理论结果是错误的。我指出了错误，给了个反例。
汉诺塔:
首先编制函数m文件，用递归的方法来做。
function hanoi(n,A,B,C)
if n==1
fprintf('%s->%s\n',A,C);
else
hanoi(n-1,A,C,B);
fprintf('%s->%s\n',A,C);
hanoi(n-1,B,A,C);
end
然后调用hanoi函数：
>> hanoi(3,'A','B','C')
A->C
A->B
C->B
A->C
B->A
B->C
A->C
