本文源於一次課題作業，部分自己寫的，部分借用了網上的demo
- 牛頓迭代法（1）
x=1:0.01:2;
y=x.^3-x.^2+sin(x)-1;
plot(x,y,'linewidth',2);grid on;%由圖像可知 根在1.05到1.15之間
syms x
s0=diff(x^3-x^2+sin(x)-1,x,1);
% 得到s0= cos(x) - 2*x + 3*x^2
% 迭代方程為 y=x-(x.^3-x.^2+sin(x)-1)/(cos(x) - 2.*x + 3*x.^2)
clear;
x=1.15;
for i=1:30
x=x-(x.^3-x.^2+sin(x)-1)/(cos(x) - 2.*x + 3*x.^2)%根據牛頓迭代法公式。一直迭代計算30次。
end
%可得 x取值為1.0935;
- 牛頓迭代法（2）
%% 繪制圖形。判斷跟的大概位置。
x=1:0.01:2;
f=x.^3-x.^2+sin(x)-1;
plot(x,f,'linewidth',2);grid on;%由圖像可知 根在1.05到1.15之間
%%
clc,clear
syms x
f=x.^3-x.^2+sin(x)-1;%所求函數
df=diff(f,x); %求取一階導數
eps=1e-4; %誤差判斷
x0=1.15; %迭代初始值。
cnt=0;
MAXCNT=20; %最大循環次數
while cnt<MAXCNT %防止無限循環
x1=x0-subs(f,x,x0)/subs(df,x,x0); %去掉這個分號,可以看到迭代過程.
if (abs(x1-x0)<eps)
break;
end
x0=x1;
cnt=cnt+1;
end
if cnt==MAXCNT
disp '不收斂'
else
vpa(x1,8)
end
- LU分解法
被調函數:
function [L,U]=lufj(A)
% 利用緊湊格式法原理 編寫的LU 分解
[n,m]=size(A); % 獲取A矩陣的行和列
if m~=n %判斷行列相等與否
error('Not a squared matrix1');
else
A(2:n)=A(2:n)/A(1,1);
for k=2:n-1
A(k,k:n)=A(k,k:n)-A(k,1:k-1)*A(1:k-1,k:n);
A(k+1:n,k)=(A(k+1:n,k)-A(k+1:n,1:k-1)*A(1:k-1,k))/A(k,k);
end %都是根據定義進行循環計算。
end
L=A;U=A;
for i=1:n
L(i,i)=1;
L(i,i+1:n)=0;
U(i,1:i-1)=0;
end
主函數：
%% 需要調用lufj函數；
A=[-2 -2 3 5
1 2 1 -2
2 5 3 -2
1 3 2 3]
b=[-1
4
7
0]
% x=A\b %左除法求解
[L,U]=lufj(A);
x0=L\b;
x=U\x0%求出的x即為解
- 拉格朗日插值法
被調函數：
function y=lagrange(x0,y0,x);
% 根據拉格朗日插值定義編寫
n=length(x0);m=length(x);
for i=1:m
z=x(i);
s=0.0;%給s的初值
for k=1:n
p=1.0;
for j=1:n
if j~=k
p=p*(z-x0(j))/(x0(k)-x0(j));
end
end
s=p*y0(k)+s;
end
y(i)=s;
end
主函數：
x=[0,1,2,4];
y=[1,9,23,3];
y0=lagrange(x,y,1.5)
- 牛頓插值
被調函數：
function yi=New_Int(x,y,xi)
%Newton基本插值公式
%x為向量，全部的插值節點
%y為向量，差值節點處的函數值
%xi為標量，是自變量
%yi為xi出的函數估計值
n=length(x);
m=length(y);
if n~=m
error('The lengths of X ang Y must be equal!');
return;
end
%計算均差表Y
Y=zeros(n);
Y(:,1)=y';
for k=1:n-1
for i=1:n-k
if abs(x(i+k)-x(i))<eps
error('the DATA is error!');
return;
end
Y(i,k+1)=(Y(i+1,k)-Y(i,k))/(x(i+k)-x(i));
end
end
%計算牛頓插值公式
yi=0;
for i=1:n
z=1;
for k=1:i-1
z=z*(xi-x(k));
end
yi=yi+Y(1,i)*z;
end
主函數：
clear all
clc
x0=[0.4 0.55 0.65 0.80 0.90 1.05];
y0=[0.41075 0.57815 0.69675 0.88811 1.0265 1.25382];
x1=0.596; % 待插值點。
y1=New_Int(x0,y0,x1)% y1即為待插值點的函數值。
TIP：主函數和被調函數要放在一個文件夾內。否則會引起調用錯誤
NOTE:本文對基本方法做了總結，你可以結合理論知識再來看代碼，希望對你有所幫助