matlab四阶龙格库塔+打靶法

matlab 专栏收录该内容
23 篇文章 1 订阅
function ys=dbf(f,a,b,a1fa,beta,h,eps)
ff=@(x,y)[y(2),f(y(1),y(2),x)];
    xvalue=a:h:b;
    n=length(xvalue)
    s0=a-0.01;
    x0=[a1fa,s0];
    flag=0;
    y0=rk4(ff,a,x0,h,a,b);
    if abs(y0(1,n)-beta)<=eps
    flag=1;
    y1=y0;
    else
        s1=s0+1;
        x0=[a1fa,s1];
        y1=rk4(ff,a,x0,h,a,b);
         if abs(y1(1,n)-beta)<=eps
             flag=1;
    end
    end
if flag~=1
    while abs(y1(1,n)-beta)>eps
        s2=s1-(y1(1,n)-beta)*(s1-s0)/(y1(1,n)-y0(1,n));
        x0=[a1fa,s2];
        y2=rk4(ff,a,x0,h,a,b);
        s0=s1;
        s1=s2;
        y0=y1;
        y1=y2;
    end
end
xvalue=a:h:b;
yvalue=y1(1,:);
ys=[xvalue',yvalue'];
function x=rk4(f,t0,x0,h,a,b)
t=a:h:b;
m=length(t);
t(1)=t0;
x(:,1)=x0;
for i=1:m-1
    L1=f(t(i),x(:,i)); 
    L2=f(t(i)+h/2,x(:,i)'+(h/2)*L1);
    L3=f(t(i)+h/2,x(:,i)'+(h/2)*L2);
      L4=f(t(i)+h,x(:,i)'+h*L3);
      x(:,i+1)=x(:,i)'+(h/6)*(L1+2*L2+2*L3+L4);
end



    
f=@(x,y,z)(x^2+z*x^2);
a=1;
b=2;
y01=0,
y0u=2,
eps=0.01;
x01=0;
x0u=2*exp(-1);
alfa=0;
beta=2;
h=0.01;
dbf(f,x01,x0u,y01,y0u,h,1e-6);
y=ans(:,2);
x=ans(:,1);
plot(x,y,'-r')
m=0:0.01:2;
n=m.*exp(-1/2*m);
plot(n,m)
plot(x,y,'-r',n,m,'-b')

龙格库塔算法是用来解微分方程的,一般学科内用于解前向的放大方程.是一个初始值问题(IVP),初始值不同解的准确度也不同.打靶法是解决边界值问题(BVP),Newton-Raphson 的shooting method 通过逐渐逼近的方法来确定初始值的算法,relaxing method 来修正以满足收敛条件和边界条件.

matlab数值计算

MATLAB中floor、round、ceil、fix区别

Matlab取整函数有: fix, floor, ceil, round.具体应用方法如下:
fix朝零方向取整,如fix(-1.3)=-1; fix(1.3)=1;
floor,顾名思义,就是地板,所以是取比它小的整数,即朝负无穷方向取整,如floor(-1.3)=-2; floor(1.3)=1;floor(-1.8)=-2,floor(1.8)=1
ceil,与floor相反,它的意思是天花板,也就是取比它大的最小整数,即朝正无穷方向取整,如ceil(-1.3)=-1; ceil(1.3)=2;ceil(-1.8)=-1,ceil(1.8)=2
round四舍五入到最近的整数,如round(-1.3)=-1;round(-1.52)=-2;round(1.3)=1;round(1.52)=2。

  • 0
    点赞
  • 4
    评论
  • 14
    收藏
  • 一键三连
    一键三连
  • 扫一扫,分享海报

相关推荐
©️2020 CSDN 皮肤主题: 技术黑板 设计师:CSDN官方博客 返回首页
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值