3point.m

三点变步长龙格库特法解非线性方程 三点变步长龙格库特法解非线性方程 三点变步长龙格库特法解非线性方程

function df=ThreePoint(func,x0,type,h)
if nargin == 3 
    h = 0.1;
else if (nargin == 4 && h == 0.0)
        disp('h不能为0！');
        return;
    end
end

y0 = subs((func), findsym(sym(func)),x0);  
y1 = subs((func), findsym(sym(func)),x0+h);
y2 = subs((func), findsym(sym(func)),x0+2*h);
y_1 = subs((func), findsym(sym(func)),x0-h);
y_2 = subs((func), findsym(sym(func)),x0-2*h);

switch type
    case 1,
        df = (-3*y0+4*y1-y2)/(2*h);     %用第一个公式求导数
       
    case 2,
        df = (3*y0-4*y_1+y_2)/(2*h);   %用第二个公式求导数
      
    case 3,
        df = (y1-y_1)/(2*h);             %用第三个公式求导数
      
end





/*变步长龙格库法求一阶微分方程*/
#include "stdio.h"
#include "math.h"

double f(double,double);
main()
{
	double y0=1,k1,k2,k3,k4,x0=0,h=0.1,h1=0.05,q,tol=0.1,y00,y01;
	int i,n=10;
	
	k1=f(x0,y0);
	k2=f(x0+h/2,y0+h/2*k1);
	k3=f(x0+h/2,y0+h/2*k2);
	k4=f(x0+h/2,y0+h*k3);
	y0=y0+h/6*(k1+2*k2+2*k3+k4);
	y00=y0;
	printf("the y00 is: %lf\n",y00);            /*步长为h时从x0出发求一步得y00*/
    k1=f(x0,y0);
	k2=f(x0+h1/2,y0+h1/2*k1);
	k3=f(x0+h1/2,y0+h1/2*k2);
	k4=f(x0+h1/2,y0+h1*k3);
	y0=y0+h1/6*(k1+2*k2+2*k3+k4);
	x0=x0+h1;
	k1=f(x0,y0);
	k2=f(x0+h1/2,y0+h1/2*k1);
	k3=f(x0+h1/2,y0+h1/2*k2);
	k4=f(x0+h1/2,y0+h1*k3);
	y0=y0+h1/6*(k1+2*k2+2*k3+k4);
	y01=y0;
	printf("the y01 is: %lf\n",y01);             /*步长为h/2时从x0出发求俩步得y01*/

	q=y01-y00;                            /*求y01和y00的差q*/
	printf("the q is: %lf\n",q);
	if(q>=tol)        /*如果q≥允许的误差tol,则反复将步长h减半，直至小于tol为止*/
	{
		while(q>=tol)
		{
			h=h/2;
            k1=f(x0,y0);
	        k2=f(x0+h/2,y0+h/2*k1);
	        k3=f(x0+h/2,y0+h/2*k2);
	        k4=f(x0+h/2,y0+h*k3);
	        y0=y0+h/6*(k1+2*k2+2*k3+k4);
	        y00=y0;
	        
	        h1=h1/2;
            k1=f(x0,y0);
	        k2=f(x0+h1/2,y0+h1/2*k1);
	        k3=f(x0+h1/2,y0+h1/2*k2);
	        k4=f(x0+h1/2,y0+h1*k3);
	        y0=y0+h1/6*(k1+2*k2+2*k3+k4);
	        x0=x0+h1;
	        k1=f(x0,y0);
	        k2=f(x0+h1/2,y0+h1/2*k1);
	        k3=f(x0+h1/2,y0+h1/2*k2);
	        k4=f(x0+h1/2,y0+h1*k3);
	        y0=y0+h1/6*(k1+2*k2+2*k3+k4);
	        y01=y0;
	        q=y01-y00;
	    }
	    
        x0=0;
	    y0=0;
	   
	    for(i=1;i<=n;i++)                /*q小于tol，用最后一次步长套用公式解题*/
	    {
	    	printf("the q, h,x0,y0 is: %lf\t%lf\t%lf\t%lf\n",q,h,x0,y0);

	    	k1=f(x0,y0);
	        k2=f(x0+h/2,y0+h/2*k1);
	        k3=f(x0+h/2,y0+h/2*k2);
	        k4=f(x0+h/2,y0+h*k3);
	      
	        y0=y0+h/6*(k1+2*k2+2*k3+k4);
            printf("%lf\t%lf\t\n",x0,y0);
            x0=x0+h;
        }
      }
      Else           /*如果q＜允许的误差tol,则反复将步长h加倍，直至大于tol为止*/
      {
      	while(q<tol)
      	{
      		h=2*h;
      		k1=f(x0,y0);
	        k2=f(x0+h/2,y0+h/2*k1);
	        k3=f(x0+h/2,y0+h/2*k2);
	        k4=f(x0+h/2,y0+h*k3);
	        y0=y0+h/6*(k1+2*k2+2*k3+k4);
            y00=y0;
            
            h1=2*h1;
            k1=f(x0,y0);
	        k2=f(x0+h1/2,y0+h1/2*k1);
	        k3=f(x0+h1/2,y0+h1/2*k2);
	        k4=f(x0+h1/2,y0+h1*k3);
	        y0=y0+h1/6*(k1+2*k2+2*k3+k4);
	        x0=x0+h1;
	        k1=f(x0,y0);
	        k2=f(x0+h1/2,y0+h1/2*k1);
	        k3=f(x0+h1/2,y0+h1/2*k2);
	        k4=f(x0+h1/2,y0+h1*k3);
	        y0=y0+h1/6*(k1+2*k2+2*k3+k4);
            y01=y0;
            q=y01-y00;
        }
        h=h/2;
        for(i=1;i<=n;i++)                /*q大于tol，用前一次一次步长套用公式解题*/
        {
        	x0=0;
        	y0=0;
        	k1=f(x0,y0);
	        k2=f(x0+h/2,y0+h/2*k1);
	        k3=f(x0+h/2,y0+h/2*k2);
	        k4=f(x0+h/2,y0+h*k3);
	        y0=y0+h/6*(k1+2*k2+2*k3+k4);
            printf("%lf\t%lf\t\n",x0,y0);
            x0=x0+h;
        }
      }
   } 
   
double f(double x,double y)                 /*定义微分方程函数*/
{
	double z;
	z=y-2*x/y;
	return z;
}