1. LU分解:
[L U]=lu(A);
2. 追赶法
function x=zhuiganfa(A,b)
[n,n]=size(A);
for i=1:n
if(i==1)
l(i)=A(i,i);
y(i)=b(i)/l(i);
else
l(i)=A(i,i)-A(i,i-1)*u(i-1);
y(i)=(b(i)-y(i-1)*A(i,i-1))/l(i);
end
if(i
u(i)=A(i,i+1)/l(i);
end
end
x(n)=y(n)
for j=n-1:-1:1
x(j)=y(j)-u(j)*x(j+1);
end
数值试验:
n=101;
关键是如何定义上述矩阵:
>>n=101;
c1=ones(1,n-1);
a1=diag(c1,-1); 这个-1说明行位置-1
c2=12*ones(1,n);
a2=diag(c2);
c3=ones(1,n-1);
a3=diag(c3,1);
a=a1+a2+a3;
3. 拉格朗日插值
function yh=lage(x,y,xh)
n=length(x);
m=length(xh);
yh=zeros(1,m);
c1=ones(n-1,1);
c2=ones(1,m);
for i=1:n
xp=x([1:i-1 i+1:n]);
yh=yh+y(i)*prod((c1*xh-xp'*c2)./(x(i)-xp'*c2));
end
end
>> x=[11,12];
>> y=[2,4];
>> xh=[11.75];
>> lage(x,y,xh)
ans =
3.5000
4 最小二乘法
1. 最小二乘的xi和yi为:
要拟合的函数为:y=a+bx-cxy 注意不是多项式
2. 编程函数为:
function z = erchen(x,y)
x1=ones(5,1);
A=[x1,x,-x.*y]; 注意点乘
z=A\y; 注意左除
a=z(1);
b=z(2);
c=z(3);
end
输入:
?? x=[1953 1964 1982 1990 2000]';
?? y=[5.82 6.95 10.08 11.34 12.66]';
?? erchen(x,y)
ans =
2.9456 =a
-0.0014 =b
-0.0005 =c
1. 最小二乘的xi和yi为:
要拟合的函数为:y=a+bx+cx2 是多项式
2. 编程函数为:
function z = erchen2(x,y)
x1=ones(5,1);
A=[x1,x,x.^2];
z=A\y;
a=z(1);
b=z(2);
c=z(3);
end
输入:
?? x=[0 0.25 0.5 0.75 1.00]';
?? y=[1.00 1.284 1.6487 2.1170 2.7183]';
?? erchen2(x,y)
ans =
1.0051
0.8642
0.8437
最小二乘多项式拟合的简单函数方法:
?? x=[0 0.25 0.5 0.75 1.00]';
?? y=[1.00 1.284 1.6487 2.1170 2.7183]';
?? P=polyfit(x,y,2) 要拟合成4次,则2改成4就可以了
P =
0.8437 0.8642 1.0051 注意此内置函数输出的结果c,b,a是反的
5 复合辛普森公式求解积分
先定义函数:
function v=f(x)
v=sin(x); “若定义有除数要点除,分母有0时要特殊定义”
end
定义程序:
function I=fsps(f,a,b,n)
h=(b-a)/n;
x=linspace(a,b,2*n+1);
y=feval(f,x);
I=(h/6)*(y(1)+2*sum(y(3:2:2*n-1))+4*sum(y(2:2:2*n))+y(2*n+1));
end
>> fsps('f ',0,1,4)
ans =
0.4597
6.不动点迭代思路
不动点迭代常常有好几个迭代的不动点函数,所以要分别定义这些函数是很困难的,如是乎使用SWITCH内置函数进行切换,叫切换函数.
1.先定义函数后进行编程的方法
先需要定义不动点函数:
function v=f(x)
v=x^3-x-1;
end
再定义编程:
function [it,x]=fixpnt1(f,a,maxit,tol)
it=0;
x=feval(f,a);
while it<=maxit & abs(x-a)>tol,
it=it+1;
a=x;
x=feval(f,a);
end
此函数的调用:
>> fixpnt1('f',2,100,1e-5)
ans =
1
3.利用切换函数SWITCH的方法(多个不动点迭代函数)
function [x,it]=fixpnt(np,a,maxit,tol)
switch np,
case 1,
phi=inline('(3*x+10)^(1/5)');
case 2,
phi=inline('sin(10*x)+2*cos(x)-3');
case 3,
phi=inline('3-atan(x)');
case 4,
phi=inline('-2-1/log(x^2+x+1)');
end
it=0;
x=phi(a);
while it<=maxit & abs(x-a)>tol,
it=it+1;
a=x;
x=phi(a);
end
使用与输入:
>> fixpnt(2,1,100,1e-5)
ans =
-4.2696
7. 雅可比迭代
function [x it]=jacobi(A,b,tol)
D=diag(diag(A));
L=D-tril(A);
U=D-triu(A);
x=zeros(size(b));
for it=1:500
x=D\(b+L*x+U*x);
error=norm(b-A*x)/norm(b);
if (error
break;
end
end
8. 高斯迭代
function [x it]=gaosi(A,b,tol)
D=diag(diag(A));
L=D-tril(A);
U=D-triu(A);
x=zeros(size(b));
for it=1:500
x=(D-L)\(b+U*x);
error=norm(b-A*x)/norm(b);
if (error
break;
end
end
9. SOR迭代
function [x it]=SOR(A,b,w,tol)
D=diag(diag(A));
L=D-tril(A);
U=D-triu(A);
x=zeros(size(b));
for it=1:500
x=(D-w*L)\(w*b+(1-w)*D*x+w*U*x);
error=norm(b-A*x)/norm(b);
if (error
break;
end
end
10.二分法:
1. 先要定义所求的函数:
function v=f(x)
v=x^3-x-1;
end
2. 二分法程序如下:
function [x,it]=erfengfa(a,b,f,tol)
fa=feval(f,a);
fb=feval(f,b);
it=0;
while abs(b-a)>tol,
it=it+1;
x=a/2+b/2;
fx=feval(f,x);
if sign(fx)==sign(fa),
a=x;fa=fx;
else
b=x;fb=fx;
end
end
11. 牛顿法:
1.先定义函数后进行编程的方法
先需要定义不动点函数
需要计算的函数f
function v = f( x )
v=x^5-3*x-10;
end
需要计算的函数的导数g
function v = g( x )
v=5*x^4-3;
end
2. 再定义编程:
function v = newton(a,f,g,maxit,tol)
it=0;
x=a;
while it<=maxit & abs(feval(f,x))>tol,
it=it+1;
x=x-feval(f,x)/feval(g,x);
end
v=[x,it];
end
12. 牛顿下山法:
1.先定义函数后进行编程的方法
先需要定义不动点函数
需要计算的函数f
function v = f(x)
v=x^2+sin(10*x)-1
end
需要计算的函数的导数g
function v =g(x)
v = 2*x+10*cos(10*x)
end
2. 再定义编程1:
function v = newtonxiashang(x0,f,g,maxit,tol)
x=x0;
it=0;
while it<=maxit & abs(feval(f,x))>tol,
it=it+1;
d=-feval(f,x)/feval(g,x);
lambda=1;
isdone=0;
while ~isdone,
xn=x+lambda*d;
if abs(feval(f,xn))
isdone=1;
else lambda=lambda*0.5;
end
end
x=xn;
end
v=[x,it];
end
3. 再定义编程2:
function v = newtonxiashang2(x0,f,g,maxit,tol)
x=x0;
it=0;
while it<=maxit & abs(feval(f,x))>tol,
it=it+1;
d=-feval(f,x)/feval(g,x);
lambda=1;
while abs(feval(f,x+lambda*d))>=abs(feval(f,x)),
lambda=0.5*lambda;
end
x=x+lambda*d;
end
v=[x,it];
end
13. 割线法
1.先定义函数后进行编程的方法
先需要定义函数
需要计算的函数f
function v = f( x )
v=x^5-3*x-10;
end
2. 再定义编程:
function v = gexian(a,b,f,maxit,tol)
it=0;
x0=a;
x=b;
while it<=maxit & abs(feval(f,x))>tol,
it=it+1;
xt=x-feval(f,x)/(feval(f,x)-feval(f,x0))*(x-x0);
x0=x;
x=xt;
end
v=[x,it];
end
输入:
>> gexian(1,2,'f',100,1e-5)
ans =
1.7226 7.0000
14. 乘幂法
function [t,y] = chenmifa( a,xinit,ep )
v0=xinit;
[tv ti]= max(abs(v0));
lam0=v0(ti);
u0=v0/lam0;
flag=0;
while (flag==0)
v1=a*u0;
[tv ti]= max(abs(v1));
lam1=v1(ti);
u0=v1/lam1;
if (abs(lam0-lam1)<=ep)
flag=1;
end
lam0=lam1;
end
t=lam1;
y=u0;
end
运行及其结果:
>> a=[2,3,2;10,3,4;3,6,1];
>> xinit=[1 1 1]';
>> ep=0.0001;
>> [t y]=chenmifa(a,xinit,ep)
t =
11.0000
y =
0.5000
1.0000
0.7500
15. 反幂法
function [t y] =fanmifa( a,xinit,ep )
v0=xinit;
[tv ti]= max(abs(v0));
lam0=v0(ti);
u0=v0/lam0;
flag=0;
while (flag==0)
v1=a^-1*u0;
[tv ti]= max(abs(v1));
lam1=v1(ti);
u0=v1/lam1;
if (abs(lam0^-1-lam1^-1)<=ep)
flag=1;
end
lam0=lam1;
end
t=lam1^-1;
y=u0;
end
运行提示:
>> a=[12 6 -6;6 16 2;-6 2 16];
>>xinit=[1 -0.5 0.5]';
>> [t y] =fanmifa( a,xinit,0.0001 )
t =
4.4560
y =
1.0000
-0.6287
0.6287
16. 结合原点平移反幂法
function [t y] =yuandian( a,p,xinit,ep )
[n n]=size(a);
v0=xinit;
[tv ti]= max(abs(v0));
lam0=v0(ti);
u0=v0/lam0;
flag=0;
while (flag==0)
v1=(a-p*eye(n))^-1*u0;
[tv ti]= max(abs(v1));
lam1=v1(ti);
u0=v1/lam1;
if (abs(lam0^-1-lam1^-1)<=ep)
flag=1;
end
lam0=lam1;
end
t=lam1^-1+p;
y=u0;
end
运行结果:
>> a=[6 2 1;2 3 1;1 1 1];
>> p=6;
>> ep=0.0001;
>> xinit=[1 1 1]';
>> [t y]=yuandian(a,p,xinit,ep)
t =
7.2880
y =
1.0000
0.5229
0.2422
17. 欧拉公式初值问题
1. 先定义导数函数f:
function v = f( x,y)
v=-y+x+1;
end
2. 定义欧拉公式编程:
function [x,y] = oula( f,y0,a,b,n )
y(1)=y0;
h=(b-a)/n;
x=a:h:b;
for i=1:n;
y(i+1)=y(i)+h*feval(f,x(i),y(i));
end
end
18. 改进的欧拉公式初值问题
1. 先定义导数函数f:
function v = f( x,y)
v=-y+x+1;
end
要先运行一下。
2. 定义欧拉公式编程:
function [x,y] =gaijingoula(f,y0,a,b,n)
y(1)=y0;
h=(b-a)/n;
x=a:h:b;
for i=1:n;
yp=y(i)+h*feval(f,x(i),y(i));
yc=y(i)+h*feval(f,x(i+1),yp);
y(i+1)=0.5*(yc+yp);
end
end
19. 梯形公式初值问题
1. 先定义导数函数f:
function v = f( x,y)
v=-y+x+1;
end
要先运行一下。
2. 定义欧拉公式编程:
function [x,y] = tixing( f,y0,a,b,n )
y(1)=y0;
h=(b-a)/n;
x=a:h:b;
for i=1:n;
y(i+1)=1/(0.5*h+1)*(y(i)+0.5*h*(feval(f,x(i),y(i))+x(i+1)+1));
end
end
20. 标准四阶四段龙格库塔公式初值问题
1. 先定义导数函数f:
function v = f( x,y)
v=y+x;
end
要先运行一下。
2. 定义欧拉公式编程:
function [x,y] = longgekuta(f,y0,a,b,n)
y(1)=y0;
h=(b-a)/n;
x=a:h:b;
for i=1:n,
k1=h*feval(f,x(i),y(i));
k2=h*feval(f,x(i)+0.5*h,y(i)+0.5*k1);
k3=h*feval(f,x(i)+0.5*h,y(i)+0.5*k2);
k4=h*feval(f,x(i)+h,y(i)+k3);
y(i+1)=y(i)+(k1+2*k2+2*k3+k4)/6;
end
end
21. 基本问题
1. 特征值:
eig(A)
2. 数值显示位数:
>> format long
>> format short
3. 左除和右除:
AX=B,则X=A\B; 左除\
XB=A,则X=A/B; 右除/ (右除是真除)
4. 特殊符号
In(x)=log(x)
根号用sqrt(x)
5.绘图功能
Plot(x,y,style)
>>x=0:0.5:10;
>>y=x;
>>hold on;
>>plot(x,y,’g+:’)
后面style的意义:
y,g,r,b分别指颜色黄色、绿色、红色、蓝色;
o * + p 代表节点的形状,圆形、*号、+号、五角星;
--、-、:代表线型,分别为虚线、线段、点线