F2
程序员文章站
2024-03-24 12:28:40
...
扫描法
% 扫描法
function p=Scanning(f,a,b,h)
c=a;i=0;
while (c<b)
f1=subs(sym(f),findsym(sym(f)),c);
f2=subs(sym(f),findsym(sym(f)),c+h);
if f1*f2<=0
i=i+1;
p(i,1)=c;
p(i,2)=c+h;
end
c=c+h;
end
对分法
%对分法
function root=HalfInterval(f,a,b,eps)
while(b-a>eps)
c=(a+b)/2;
f1=subs(sym(f),findsym(sym(f)),c);
f2=subs(sym(f),findsym(sym(f)),a);
if f1==0
break;
elseif f1*f2<0
b=c;
else
a=c;
end
end
root=c;
立方根
%立方根
function root=CubeRoot(a,eps)
to1=1;x=a;
root=x-(x^3-a)/(3*x^2);
while(to1>eps)
x=root;
root=x-(x^3-a)/(3*x^2);
to1=abs(root-x);
end
牛顿法
%牛顿法
function root=NEWtonRoot(f,x,eps)
to1=1;
fun=diff(sym(f));
fa=subs(sym(f),findsym(sym(f)),x);
dfa=subs(sym(fun),findsym(sym(fun)),x);
root=x-fa/dfa;
while (to1>eps)
x=root;
fa=subs(sym(f),findsym(sym(f)),x);
dfa=subs(sym(fun),findsym(sym(fun)),x);
root=x-fa/dfa;
to1=abs(root-x);
end
弦割法
%弦割法
function root=Secant(f,a,b,eps)
tol=1;
fa=subs(sym(f),findsym(sym(f)),a);
fb=subs(sym(f),findsym(sym(f)),b);
root=b-(b-a)*fb/(fb-fa);
while (tol>eps)
a=b;
b=root;
fa=subs(sym(f),findsym(sym(f)),a);
fb=subs(sym(f),findsym(sym(f)),b);
root=b-(b-a)*fb/(fb-fa);
tol=abs(root-b);
end
高斯列主元
%高斯列主元
function x=GaussMain(A,b)
N=size(A);
n=N(1);
index=0;
A(:,n+1)=b;
for i=1:(n-1)
me=max(abs(A(i:n,i)));
for k=i:n
if(abs(A(k,i))==me)
index=k;
break
end
end
if i~=index
temp=A(i,i:n+1);
A(i,i:n+1)=A(index,i:n+1);
A(index,i:n+1)=temp;
end
for j=(i+1):n
L=A(j,i);
m=A(i,i);
A(j,i:n+1)=A(j,i:n+1)-L*A(i,i:n+1)/m;
end
end
x=SolveUpTriangle(A);
function x=SolveUpTriangle(A)
N=size(A);
n=N(1);
x(n)=A(n,n+1)/A(n,n);
for i=n-1:-1:1
x(i)=(A(i,n+1)-sum(A(i,i+1:n).*x(i+1:n)))/A(i,i);
end
高斯约当
%高斯约当
function x=GaussJordan(A,b)
N=size(A);
n=N(1);
A(:,n+1)=b;
index=0;
for i=1:1:n
me=max(abs(A(i:n,i)));
for k=i:n
if(abs(A(k,i))==me)
index=k;
break
end
end
if i~=index
temp=A(i,i:n+1);
A(i,i:n+1)=A(index,i:n+1);
A(index,i:n+1)=temp;
end
A(i,i:n+1)=A(i,i:n+1)/A(i,i);
for j=1:1:n
if j~=i
A(j,i:n+1)=A(j,i:n+1)-A(j,i)*A(i,i:n+1);
end
end
end
x=A(1:n,n+1);
LU分解
%LU分解
function x=Doolittle(A,b)
N=size(A);
n=N(1);
L=eye(n,n);
U=zeros(n,n);
U(1,1:n)=A(1,1:n);
L(1:n,1)=A(1:n,1)/U(1,1);
for k=2:n
for i=k:n
U(k,i)=A(k,i)-L(k,1:(k-1))*U(1:(k-1),i);
end
for j=(k+1):n
L(j,k)=(A(j,k)-L(j,1:(k-1))*U(1:(k-1),k))/U(k,k);
end
end
y=SolveDownTriangle(L,b);
x=SolveUpTriangle(U,y);
function x=SolveDownTriangle(A,b)
N=size(A);
n=N(1);
x(1)=b(1);
for i=2:1:n
x(i)=b(i)-sum(A(i,1:i-1).*x(1:i-1));
end
function x=SolveUpTriangle(A,b)
N=size(A);
n=N(1);
x(n)=b(n)/A(n,n);
for i=n-1:-1:1
x(i)=(b(i)-sum(A(i,i+1:n).*x(i+1:n)))/A(i,i);
end
追赶法
function x=followup(A,d)
n=rank(A);
b=ones(n,1);
a=ones(n,1);
c=ones(n-1,1);
for i=1:n-1
a(i+1)=A(i+1,i);
c(i)=A(i,i+1);
b(i)=A(i,i);
end
b(n)=A(n,n);
c(1)=c(1)/b(1);
d(1)=d(1)/b(1);
for k=2:1:n-1
c(k)=c(k)/(b(k)-a(k)*c(k-1));
d(k)=(d(k)-a(k)*d(k-1))/(b(k)-a(k)*c(k-1));
end
d(n)=(d(n)-a(n)*d(n-1))/(b(n)-a(n)*c(n-1));
x(n)=d(n);
for k=n-1:-1:1
x(k)=d(k)-c(k)*x(k+1);
end
Jacobian迭代
%Jacobian迭代
function [x,n]=Jacobi(A,b,x0,eps,M)
D=diag(diag(A));
L=-tril(A,-1);
U=-triu(A,1);
B=D\(L+U);
f=D\b;
x=B*x0+f;
n=1;
while norm(x-x0)>eps
x0=x;
x=B*x0+f;
n=n+1;
if(n>=M)
disp('Warning: 迭代式未收敛');
return;
end
end
Gauss-Seridel迭代
%Gauss-Seridel迭代
function [x,n]=Gauseidel(A,b,x0,eps,M)
D=diag(diag(A));
L=-tril(A,-1);
U=-triu(A,1);
G=(D-L)\U;
f=(D-L)\b;
x=G*x0+f;
n=1;
while norm(x-x0)>=eps
x0=x;
x=G*x0+f;
n=n+1;
if(n>=M)
disp('Warning :');
return;
end
end
乘幂法
%乘幂法
function [l,v,s]=pmethod(A,x0,eps)
v=x0;
M=5000;
m=0;
l=0;
for k=1:1:M
y=A*v;
m=absmaxcomponent(y);
v=y/m;
if abs(m-l)<eps
l=m;
s=k;
return;
else
if k==M
disp('迭代步数太多，收敛速度太慢！');
l=m;
s=M;
else
l=m;
end
end
end
function me=absmaxcomponent(y)
n=length(y);
index=0;
me=max(abs(y));
for k=1:1:n
if(abs(y(k))==me)
index=k;
break;
end
end
me=y(index);
反幂法
%反幂法
function [l,v,s]=ipmethod(A,x0,eps)
v=x0;
M=5000;
m=0;
l=0;
for k=1:1:M
y=A\v;
m=absmaxcomponent(y);
v=y/m;
if abs(m-l)<eps
l=1/m;s=k;return;
else
if k==M
disp('迭代步数太多，收敛速度太慢！');
l=1/m;
s=M;
else
l=m;
end
end
end
function me=absmaxcomponent(y)
n=length(y);
index=0;
me=max(abs(y));
for k=1:1:n
if(abs(y(k))==me)
index=k;
break;
end
end
me=y(index);
Lagrange插值
%Lagrange插值
function f0 = lagrange(x,y,x0)
f0=0.0;
n=length(x);
for i=1:n
L=y(i);
for j=1:n
if j~=i
L=L*(x0-x(j))/(x(i)-x(j));
end
end
f0=f0+L;
end
Netwon插值
%Netwon插值
function f0=Newton(x,y,x0)
n=length(x);
for k=1:1:n-1
for i=n:-1:k+1
y(i)=(y(i)-y(i-1))/(x(i)-x(i-k));
end
end
f0=y(1);
h=1;
for i=2:1:n
h=h*(x0-x(i-1));
f0=f0+h*y(i);
end
Newton前向插值
%Netwon前向
function f0 = Newtonforward(x,y,x0)
n=length(x);
c(1:n)=0.0;
f0=y(1);
y1=0;
for i=1:n-1
for j=1:n-i
y1(j)=y(j+1)-y(j);
end
c(i)=y1(1);
L=(x0-x(1))/(x(2)-x(1));
for k=1:i-1
L=L*((x0-x(1))/(x(2)-x(1))-k);
end
f0=f0+c(i)*L/factorial(i);
y=y1;
end
Newton后向插值
%Newton后向插值
function f0=Nwetonback(x,y,x0)
n=length(x);
c(1:n)=0.0;
f0=y(n);
y1=0;
for i=1:n-1
for j=i+1:n
y1(j)=y(j)-y(j-1);
end
c(i)=y1(n);
L=(x0-x(n))/(x(2)-x(1));
for k=1:i-1
L=L*((x0-x(n))/(x(2)-x(1))+k);
end
f0=f0+c(i)*L/factorial(i);
y=y1;
end
分段Hermite插值
%分段Hermite插值
function f0=SubHermite(x,y,y_1,x0)
n=length(x);
for i=1:n
if(x(i)<=x0)&&(x(i+1)>=x0)
index=i;
break;
end
end
h=x(index+1)-x(index);
f1=y(index)*(1+2*(x0-x(index))/h)*(x0-x(index+1))^2/h/h+...
y(index+1)*(1-2*(x0-x(index+1))/h)*(x0-x(index))^2/h/h;
fr=y_1(index)*(x0-x(index))*(x0-x(index+1))^2/h/h+...
y_1(index+1)*(x0-x(index+1))*(x0-x(index))^2/h/h;
f0=f1+fr;
曲线拟合
%曲线拟合
function A=multifit(x,y,m)
N=length(x);
M=length(y);
c(1:(2*m+1))=0;
b(1:(m+1))=0;
for j=1:(2*m+1)
for k=1:N
c(j)=c(j)+x(k)^(j-1);
if(j<(m+2))
b(j)=b(j)+y(k)*x(k)^(j-1);
end
end
end
C(1,:)=c(1:(m+1));
for s=2:(m+1)
C(s,:)=c(s:(m+s));
end
A=inv(C)*b';
复合梯形
%复合梯形
function I=FHTrapr1(f,a,b,N)
%f 被积函数
%a，b 积分上下限
%N 区间等分数
%I 积分结果
h=(b-a)/N;
s=0;
for i=1:N-1
x=a+i*h;
s=s+subs(sym(f),findsym(sym(f)),x);
end
I=(h/2)*(subs(sym(f),findsym(sym(f)),a)+...
2*s+subs(sym(f),findsym(sym(f)),b));
复合Simpson
%复合Simpson
function I=FHSimpson(f,a,b,N)
%f 被积函数
%a,b 积分上下限
%N 区间等分数
%I 积分结果
h=(b-a)/N;
s1=0;
s2=0;
for i=1:N
x=a+(2*i-1)*(h/2);
s1=s1+subs(sym(f),findsym(sym(f)),x);
end
for i=1:N-1
x=a+(2*i)*(h/2);
s2=s2+subs(sym(f),findsym(sym(f)),x);
end
I=(h/6)*(subs(sym(f),findsym(sym(f)),a)+...
4*s1+2*s2+subs(sym(f),findsym(sym(f)),b));
复合Cotes
%复合Cotes
function I=FHCotes(f,a,b,N)
%f 被积函数
%a,b 积分上下限
%N 区间等分数
%I 积分结果
h=(b-a)/N;
s1=0;s2=0;s3=0;s4=0;
for i=1:N
x=a+(4*i-3)*(h/4);
s1=s1+subs(sym(f),findsym(sym(f)),x);
end
for i=1:N
x=a+(4*i-2)*(h/4);
s2=s2+subs(sym(f),findsym(sym(f)),x);
end
for i=1:N
x=a+(4*i-1)*(h/4);
s3=s3+subs(sym(f),findsym(sym(f)),x);
end
for i=1:N-1
x=a+(4*i)*(h/4);
s4=s4+subs(sym(f),findsym(sym(f)),x);
end
I=(h/90)*(7*subs(sym(f),findsym(sym(f)),a)...
+32*s1+12*s2+32*s3+14*s4+7*subs(sym(f),findsym(sym(f)),b));
变步长梯形
%变步长梯形
function [I,step] = ChangeTrapr1(f,a,b,eps)
%f 被积函数
%a,b 积分上下限
%eps 精度
%I 积分结果
%step 积分的子区间数
if(nargin==3)
eps=1.0e-4;
end
k=0;
h=b-a;
T1=0;
T2=(h/2)*(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b));
while abs(T2-T1)>eps
T1=T2;
h=h/2;
Q=0;
for i=1:2^k
x=a+h*(2*i-1);
Q=Q+subs(sym(f),findsym(sym(f)),x);
end
T2=T1/2+h*Q;
k=k+1;
end
I=T2;
step=2^k;
Romberg积分
%Romberg积分
function [I,step]=Roberg(f,a,b,eps)
if(nargin==3)
eps=1.0e-4;
end
M=1;
tol=10;
k=0;
T=zeros(1,1);
h=b-a;
T(1,1)=(h/2)*(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b));
while tol>eps
k=k+1;
h=h/2;
Q=0;
for i=1:M
x=a+h*(2*i-1);
Q=Q+subs(sym(f),findsym(sym(f)),x);
end
T(k+1,1)=T(k,1)/2+h*Q;
M=2*M;
for j=1:k
T(k+1,j+1)=T(k+1,j)+(T(k+1,j)-T(k,j))/(4^j-1);
end
tol=abs(T(k+1,j+1)-T(k,j));
end
I=T(k+1,k+1);
step=2^k;
上一篇： Python：OS模块
下一篇： 图片缩放