1
创建一个函数
%牛顿法求立方根
function
x=cube_newton(a)
f=@(x)x^3-a;
df=diff(sym('x^3-a'));
if
a==0;
x1=a;
else
x0=a;
x1=x0-f(x0)/subs(df,x0);
while
abs(x1-x0)>1e-6
x0=x1;
x1=x0-f(x0)/subs(df,x0);
end
end
x=x1;
2
调用求解
>>
a=cube_newton(5)
a
=
17100
>>
%% set para
d=6;
tol=1e-5;
maxIter=100;
r=-2:001:2; %实部虚部的范围
[x y]=meshgrid(r); %产生实部虚部二维网格
Z=x+1iy; %Z对应网格的虚平面
%% Define fuction
f=@(x,d) (x^d)-1;
fprime=@(x,d) d(x^(d-1));
%% Perform Newton iterations
for k=1:maxIter;
Z=Z-(f(Z,d)/fprime(Z,d));
end
%% Find d roots of unity, and the mask
renderMat=0;
for j=1:d
root=exp(2pi1i/d)^j; % the jth root
Mj=abs(Z-root); % distance Z中每点都这个根的距离
% Each root gets a unique number in [1,d]
mask=(Mj<=tol)j; %Mj<=tol返回满足误差的逻辑矩阵
%满足误差部分为1j,不满足部分为0
renderMat=renderMat+mask;
%加起来之后renderMat中收敛于第j个根的区域数据都是j
%那么收敛于第j个根的区域都是同一种颜色
end
colormap(hsv(d+1)); % Set the color map
imagesc(r,r,renderMat) % Render the fractal
xlabel('Re(Z)');ylabel('Im(Z)');
h=colorbar;
set(h,'ytick',(2(0:d)+1)d/(d+1)/2);
str=arrayfun(@(x)num2str(x,'%2f'),exp(2pi1i/d)^(1:d),'uniformoutput',false);
set(h,'yticklabel',[{'未收敛'},str]);
在Matlab命令行下输入:
edit
Gauss_seidelm
然后将下面输入,并保存
function
x=Gauss_Seidel(A,b,x0,tol)
if
(nargin==2)
x0=ones(size(b));
tol=1e-6;
elseif
(nargin==3)
tol=1e-6;
else
sprintf('USAGE:Gauss_Seidel(A,b,x0,tol)')
end
D=diag(diag(A));
U=triu(A,1);
L=tril(A,-1);
G=-(D+L)\U;
d1=(D+L)\b;
x=Gx0+d1;
n=1;
while
norm(x-x0)>=tol
x0=x;
x=Gx0+d1;
n=n+1;
end
n
到Matlab主窗口,输入:
A=[10,-1,2,0;-1,11,-1,3;2,-1,10,-1;0,3,-1,8];
b=[6;25;-11;15];
x=gauss_seidel(A,b)
就可以得到结果。
x
=
10000
20000
-10000
10000
迭代8次。
拉格朗日function y=lagrange(x0,y0,x)
n=length(x0);m=length(x);
for i=1:m
z=x(i);
s=00;
for k=1:n
p=10;
for j=1:n
if j~=k
p=p(z-x0(j))/(x0(k)-x0(j));
end
end
s=py0(k)+s;
end
y(i)=s;
end SOR迭代法的Matlab程序
function [x]=SOR_iterative(A,b)
% 用SOR迭代求解线性方程组,矩阵A是方阵
x0=zeros(1,length(b)); % 赋初值
tol=10^(-2); % 给定误差界
N=1000; % 给定最大迭代次数
[n,n]=size(A); % 确定矩阵A的阶
w=1; % 给定松弛因子
k=1;
% 迭代过程
while k<=N
x(1)=(b(1)-A(1,2:n)x0(2:n)')/A(1,1);
for i=2:n
x(i)=(1-w)x0(i)+w(b(i)-A(i,1:i-1)x(1:i-1)'-A(i,i+1:n)x0(i+1:n)')/A(i,i);
end
if max(abs(x-x0))<=tol
fid = fopen('SOR_iter_resulttxt', 'wt');
fprintf(fid,'\n用SOR迭代求解线性方程组的输出结果\n\n');
fprintf(fid,'迭代次数: %d次\n\n',k);
fprintf(fid,'x的值\n\n');
fprintf(fid, '%128f \n', x);
break;
end
k=k+1;
x0=x;
end
if k==N+1
fid = fopen('SOR_iter_resulttxt', 'wt');
fprintf(fid,'\n用SOR迭代求解线性方程组的输出结果\n\n');
fprintf(fid,'迭代次数: %d次\n\n',k);
fprintf(fid,'超过最大迭代次数,求解失败!');
fclose(fid);
end Matlab中龙格-库塔(Runge-Kutta)方法原理及实现龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高,采取措施对误差进行抑制,所以其实现原理也较复杂。该算法是构建在数学支持的基础之上的。龙格库塔方法的理论基础来源于泰勒公式和使用斜率近似表达微分,它在积分区间多预计算出几个点的斜率,然后进行加权平均,用做下一点的依据,从而构造出了精度更高的数值积分计算方法。如果预先求两个点的斜率就是二阶龙格库塔法,如果预先取四个点就是四阶龙格库塔法。一阶常微分方程可以写作:y'=f(x,y),使用差分概念。
(Yn+1-Yn)/h= f(Xn,Yn)推出(近似等于,极限为Yn')
Yn+1=Yn+hf(Xn,Yn)
另外根据微分中值定理,存在0<t<1,使得
Yn+1=Yn+hf(Xn+th,Y(Xn+th))
这里K=f(Xn+th,Y(Xn+th))称为平均斜率,龙格库塔方法就是求得K的一种算法。
利用这样的原理,经过复杂的数学推导(过于繁琐省略),可以得出截断误差为O(h^5)的四阶龙格库塔公式:
K1=f(Xn,Yn);
K2=f(Xn+h/2,Yn+(h/2)K1);
K3=f(Xn+h/2,Yn+(h/2)K2);
K4=f(Xn+h,Yn+hK3);
Yn+1=Yn+h(K1+2K2+2K3+K4)(1/6);
所以,为了更好更准确地把握时间关系,应自己在理解龙格库塔原理的基础上,编写定步长的龙格库塔函数,经过学习其原理,已经完成了一维的龙格库塔函数。
仔细思考之后,发现其实如果是需要解多个微分方程组,可以想象成多个微分方程并行进行求解,时间,步长都是共同的,首先把预定的初始值给每个微分方程的第一步,然后每走一步,对多个微分方程共同求解。想通之后发现,整个过程其实很直观,只是不停的逼近计算罢了。编写的定步长的龙格库塔计算函数:
function [x,y]=runge_kutta1(ufunc,y0,h,a,b)%参数表顺序依次是微分方程组的函数名称,初始值向量,步长,时间起点,时间终点(参数形式参考了ode45函数)
n=floor((b-a)/h);%求步数
x(1)=a;%时间起点
y(:,1)=y0;%赋初值,可以是向量,但是要注意维数
for ii=1:n
x(ii+1)=x(ii)+h;
k1=ufunc(x(ii),y(:,ii));
k2=ufunc(x(ii)+h/2,y(:,ii)+hk1/2);
k3=ufunc(x(ii)+h/2,y(:,ii)+hk2/2);
k4=ufunc(x(ii)+h,y(:,ii)+hk3);
y(:,ii+1)=y(:,ii)+h(k1+2k2+2k3+k4)/6;
%按照龙格库塔方法进行数值求解
end
调用的子函数以及其调用语句:
function dy=test_fun(x,y)
dy = zeros(3,1);%初始化列向量
dy(1) = y(2) y(3);
dy(2) = -y(1) + y(3);
dy(3) = -051 y(1) y(2);
对该微分方程组用ode45和自编的龙格库塔函数进行比较,调用如下:
[T,F] = ode45(@test_fun,[0 15],[1 1 3]);
subplot(121)
plot(T,F)%Matlab自带的ode45函数效果
title('ode45函数效果')
[T1,F1]=runge_kutta1(@test_fun,[1 1 3],025,0,15);%测试时改变test_fun的函数维数,别忘记改变初始值的维数
subplot(122)
plot(T1,F1)%自编的龙格库塔函数效果
title('自编的 龙格库塔函数')
以上就是关于请问怎么用matlab编程,使用牛顿迭代法求根号5的立方的近似值全部的内容,包括:请问怎么用matlab编程,使用牛顿迭代法求根号5的立方的近似值、一段画牛顿迭代收敛域的matlab程序,求解释、求!!在MATLAB中用求解高斯—塞德尔迭代法线性方程组程序等相关内容解答,如果想了解更多相关内容,可以关注我们,你们的支持是我们更新的动力!
欢迎分享,转载请注明来源:内存溢出
评论列表(0条)