如何用MATLAB编写的拉格朗日插值算法的程序、二阶龙格-库塔方法的程序和SOR迭代法的程序

如何用MATLAB编写的拉格朗日插值算法的程序、二阶龙格-库塔方法的程序和SOR迭代法的程序,第1张

拉格朗日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('自编的 龙格库塔函数')

x0=-200;x1=000;x2=100;x3=200;

y0=1700;y1=100;y2=200;y3=1700;

x=06

y=(x-x1)(x-x2)(x-x3)/((x0-x1)(x0-x2)(x0-x3))y0+(x-x0)(x-x2)(x-x3)/((x1-x0)(x1-x2)(x1-x3))y1+(x-x0)(x-x1)(x-x3)/((x2-x0)(x2-x1)(x2-x3))y2+(x-x0)(x-x1)(x-x2)/((x3-x0)(x3-x1)(x3-x2))y3;

disp('y=');disp(y);

结果为:x =

06000

y=

02560

用拉格朗日插值做y=1/(1+x^2);本来就有很大的偏差,你的插值函数没写错

这是我做过的一个y=1/(1+25x^2)的图像 ,这里取得是11点,你的应该是5点的吧。你的做完是这种么。这样没错,插值法不适合太多点,想做多点的用的是分段插值。你应该是在学计算方法的课程吧,后边就会讲的。

在MATLAB中,一维多项式插值的方法通过命令interp1实现,其具体的调用格式如下:

插值的方法method参数的取值和对应的含义如下:

nearest:最邻近插值方法(nearest neighbor interpolation)。这种插值方法在已知数据的最邻近点设置插值点,对插值点的数值进行四舍五入,对超出范围的数据点返回NaN。

linear:线性插值(linear interpolation),这是interp1命令中method的默认数值。该方法采用直线将相邻的数据点相连,对超出数据范围的数据点返回NaN。

spline:三次样条插值(cubic spline interpolation),该方法采用三次样条函数获取插值数据点,在已知点为端点的情况下,插值函数至少具有相同的一阶和二阶导数。

pchip:分段三次厄米多项式差值(piecewise cubic Hermite interpolation)。

cubic:三次多项式插值,与分段三次厄米多项式插值方法相同。

v5cubic:MATLAB5中使用的三次多项式插值。

spline函数可以实现三次样条插值

x=0:10;

y=sin(x);

xx=0:25:10;

yy=spline(x,y,xx);

plot(x,y,'o',xx,yy)

另外fnpltcsapi这两个函数也是三次样条插值函数,具体你可以help一下!

现在电脑上没有matlab,一会给你程序,呵呵!

以上就是关于如何用MATLAB编写的拉格朗日插值算法的程序、二阶龙格-库塔方法的程序和SOR迭代法的程序全部的内容,包括:如何用MATLAB编写的拉格朗日插值算法的程序、二阶龙格-库塔方法的程序和SOR迭代法的程序、急用matlab编写拉格朗日插值算法的程序、matlab拉格朗日插值程序等相关内容解答,如果想了解更多相关内容,可以关注我们,你们的支持是我们更新的动力!

欢迎分享,转载请注明来源:内存溢出

原文地址:https://54852.com/zz/9511401.html

(0)
打赏 微信扫一扫微信扫一扫 支付宝扫一扫支付宝扫一扫
上一篇 2023-04-29
下一篇2023-04-29

发表评论

登录后才能评论

评论列表(0条)

    保存