四阶龙格库塔matlab实现遇到的问题

function [X,Y]=odeRK4(str,a,b,ya,n)
% 四阶龙格-库塔方法
%f 函数
%a和b 求解区域的端点
%ya 初始条件y(a)
%n 求解步数
%f 如果是文件,调用[x,y]=odeRK4(@f,a,b,ya,n).
%f 如果是匿名函数,调用[x,y]=odeRK4(f,a,b,ya,n).
h=(b-a)/n;
Y=zeros(1,n+1);
X=a:h:b;
Y(1)=ya;
f_xy= @(x,y)str;
for j=1:n
k1=f_xy(X(j),Y(j));
k2=f_xy(X(j)+h/2,Y(j)+h/2*k1);
k3=f_xy(X(j)+h/2,Y(j)+h/2*k2);
k4=f_xy(X(j)+h,Y(j)+h*k3);
Y(j+1)=Y(j)+h/6*(k1+k4)+h/3*(k2+k3);
end
for k=1:n+1
fprintf('x[%d]=%f\ty[%d]=%f\n',k-1,X(k),k-1,Y(k));
end
end

clear all
figure('name','四阶龙格-库塔', 'NumberTitle','off') %建立名字是四阶龙格-库塔的图形窗口
uicontrol(gcf,'style','text', 'position',[50 350 350 50], 'string','常微分方程组初值问题数值解', 'fontsize',16);
% 建立名为常微分方程组初值问题数值解的求解的静态编辑框
text_fun1= uicontrol(gcf, 'Style','Text','String','常微分方程表达式: ','FontSize',14,'Position',[55 278 160 25],'HorizontalAlignment','Left');
edit_fun1= uicontrol(gcf, 'Style','Edit','String','-x*x*y*y*y','Position',[220 278 140 25],'HorizontalAlignment','left', 'FontSize',14);
text_fun2= uicontrol(gcf, 'Style','Text','String','区间左端点: ','FontSize',14,'Position',[55 248 160 25],'HorizontalAlignment','Left');
edit_fun2= uicontrol(gcf,'style','Edit','String','0', 'position',[220 248 140 25],'HorizontalAlignment','left', 'FontSize',14);
text_fun3= uicontrol(gcf, 'Style','Text','String','区间右端点: ','FontSize',14,'Position',[55 218 160 25],'HorizontalAlignment','Left');
edit_fun3= uicontrol(gcf,'style','Edit','String','5', 'position',[220 218 140 25],'HorizontalAlignment','left', 'FontSize',14);
text_fun4= uicontrol(gcf, 'Style','Text','String','初值: ','FontSize',14,'Position',[55 188 160 25],'HorizontalAlignment','Left');
edit_fun4= uicontrol(gcf,'style','Edit','String','1', 'position',[220 188 140 25],'HorizontalAlignment','left', 'FontSize',14);
text_fun5= uicontrol(gcf, 'Style','Text','String','求解步数: ','FontSize',14,'Position',[55 158 160 25],'HorizontalAlignment','Left');
edit_fun5= uicontrol(gcf,'style','Edit','String','20', 'position',[220 158 140 25],'HorizontalAlignment','left', 'FontSize',14);
uicontrol(gcf,'style','push', 'position',[55 108 150 30], 'string','开始运行', 'fontsize',10, 'call',['odefun =num2str(get(edit_fun1,''string''));a = str2num(get(edit_fun2,''string'')),b = str2num(get(edit_fun3,''string'')),ya = str2num(get(edit_fun4,''string'')),n = str2num(get(edit_fun5,''string''));','[x,y]=odeRK4(odefun,a,b,ya,n)'])
% 建立名为开始运行的静态编辑框
push_Quit = uicontrol(gcf, 'Style','Push','String','Quit','value',0,'FontSize',10, 'Pos',[450 50 100 30], 'Call', 'fun_quit');

运行报错
在赋值 A(:) = B 中,A 和 B 中的元素数目必须相同。

出错 odeRK4 (line 19)
Y(j+1)=Y(j)+h/6*(k1+k4)+h/3*(k2+k3);

计算 UIControl Callback 时出错

1个回答

Csdn user default icon
上传中...
上传图片
插入图片
抄袭、复制答案,以达到刷声望分或其他目的的行为,在CSDN问答是严格禁止的,一经发现立刻封号。是时候展现真正的技术了!
其他相关推荐
求大佬帮我看一下我的四阶龙格库塔法求解常微分方程组的MATLAB程序对吗

![常微分方程组](https://img-ask.csdn.net/upload/201905/15/1557892716_658772.png) 1. ``` function f=myfun(t,y,beta1,beta2,epsilong1,epsilong2,delta,u1,u2,k,alpha,p) f(1)=-beta1.*epsilong1.*y(1).*(1-u1).*y(3)-beta2.*epsilong2.*y(1).*(1-u1).*delta.*y(2); f(2)=beta1.*epsilong1.*y(1).*(1-u1).*y(3)-beta2.*epsilong2.*y(1).*(1-u1).*delta.*y(2)-k*y(2); f(3)=k.*y(2)-alpha.*y(3)-u2.*y(3); f(4)=p.*alpha.*y(3)+u2.*y(3); f=f( : ); ``` 这是我的函数m文件 2. ``` function [ t, y]=RK4(myfun,tspan,y0,h) t=tspan(1):h:tspan(2); y=zeros(length(y0),length(t)); y(:,1)=y0(:); for n=1:(length(t)-1) k1=feval(myfun,t(n),y(:,n)); k2=feval(myfun,t(n)+h/2,y(:,n)+h/2*k1); k3=feval(myfun,t(n)+h/2,y(:,n)+h/2*k2); k4=feval(myfun,t(n+1),y(:,n)+h*k3); y(:,n+1)=y(:,n)+(h/6).*(k1+2*k2+2*k3+k4); end ``` 这是四阶龙格库塔的程序文件 3. ``` clear; clc; tspan=[0 300]; y0=[1000000 0 1 0]; % 初值 beta1=0.8; beta2=0.5; epsilong1=0.5; epsilong2=0.2; delta=0.5; u1=0.7; u2=0.03; k=0.526; alpha=0.244; p=0.93; [t, y]=RK4(@(t,y)myfun(t,y,beta1,beta2,epsilong1,epsilong2,delta,u1,u2,k,alpha,p),tspan,y0,5); plot(t,y); title('系统方程');xlabel('时间(t)');ylabel('y'); legend('y_1','y_2','y_3','y_4',2); ``` 这是最后的程序,上面附了一些初值,最后得出的结果是这样的, ![图片说明](https://img-ask.csdn.net/upload/201905/15/1557893185_84919.jpg) 这是线性的吗,感觉很懵逼啊,而且y的值怎么会出现负的

恳请指点:四阶龙格库塔解常微分方程不断改变步长,运行有时候出结果,有时候程序没有结果

d更改很多次,查了不少资料,仍没有结果,恳请大神能帮忙指点一下,或者该参考什么方向的书籍?非常感 谢,当C[0]=200;时 //功能:方程组1和方程组2解不断逼近,t不断改变方程组1,A不断改变方程组2迭代初值,误差满足要求, 输出此时A、t。现在在步长为0.01,C[0]=100;时运行出t=1.1,A=43.0,改变C[0]=200就没有结果了1、参看 资料变了很多步长,最后面的是看的相关论文改的根据误差选取不同的步长,不出结果2,Db[i]和Dc[i]同时 满足,A[i]为负值(理论上为正) #include<stdio.h> #include <math.h> #include<stdlib.h> #define PI 3.141593 #define Vx y[1] #define Vy y[2] #define X y[3] #define Y y[4] #define f1(y1,y2,y3,y4,t) (-0.1109*y1) #define f2(y1,y2,y3,y4,t) (-0.1109*y2-9.81)//方程2 #define f3(y1,y2,y3,y4,t) (y1) #define f4(y1,y2,y3,y4,t) (y2) void main(void) { int i,j,V,m=4,n=151; double t,A[150],Db[150],Dc[150],B[150],C[150],G[150],M[150],r[150],R[150],d,y[150],y0 [150],k1[150],k2[150],k3[150],k4[150]; B[0]=0; C[0]=100; G[0]=100; M[0]=0; M[1]=0; M[2]=-10; V=120; d=0.01; t=0; B[1]=B[0]+M[0]*t;//方程 C[1]=C[0]+M[1]*t; G[1]=G[0]+M[2]*t; r[1]=sqrt(B[1]*B[1]+G[1]*G[1]); R[1]=sqrt(r[1]*r[1]+C[1]*C[1]); A[1]=atan(C[1]/r[1]); A[1]=A[1]*180/PI;//弧度化角度 A[1]=cos(A[1]*PI/180);//输入必须为弧度 y0[1]=V*A[1];//初始条件 A[1]=atan(C[1]/r[1]); A[1]=A[1]*180/PI;//弧度化为度 A[1]=sin(A[1]*PI/180);//输入必须为弧度 y0[2]=V*A[1]; y0[3]=0; y0[4]=0; for(j=1;j<=m;j++) y[j]=y0[j]; for(i=1;i<n;i++) { for(j=1;j<=m;j++) { if(j==1) k1[j]=d*f1((y[j]),(y[j+1]),(y[j+2]),(y[j+3]),t); else if(j==2) k1[j]=d*f2((y[j-1]),(y[j]),(y[j+1]),(y[j+2]),t); else if(j==3) k1[j]=d*f3((y[j-2]),(y[j-1]),(y[j]),(y[j+1]),t); else if(j==4) k1[j]=d*f4((y[j-3]),(y[j-2]),(y[j-1]),(y[j]),t); } for(j=1;j<=m;j++) { if(j==1) k2[j]=d*f1((y[j]+0.5*k1[j]),(y[j+1]+0.5*k1[j+1]),(y[j+2]+0.5*k1[j+2]),(y[j +3]+0.5*k1[j+3]),(t+0.5*d)); else if(j==2) k2[j]=d*f2((y[j-1]+0.5*k1[j-1]),(y[j]+0.5*k1[j]),(y[j+1]+0.5*k1[j+1]),(y[j +2]+0.5*k1[j+2]),(t+0.5*d)); else if(j==3) k2[j]=d*f3((y[j-2]+0.5*k1[j-2]),(y[j-1]+0.5*k1[j-1]),(y[j]+0.5*k1[j]),(y[j +1]+0.5*k1[j+1]),(t+0.5*d)); else if(j==4) k2[j]=d*f4((y[j-3]+0.5*k1[j-3]),(y[j-2]+0.5*k1[j-2]),(y[j-1]+0.5*k1[j-1]),(y [j]+0.5*k1[j]),(t+0.5*d)); } for(j=1;j<=m;j++) { if(j==1) k3[j]=d*f1((y[j]+0.5*k2[j]),(y[j+1]+0.5*k2[j+1]),(y[j+2]+0.5*k2[j+2]),(y[j +3]+0.5*k2[j+3]),(t+0.5*d)); else if(j==2) k3[j]=d*f2((y[j-1]+0.5*k2[j-1]),(y[j]+0.5*k2[j]),(y[j+1]+0.5*k2[j+1]),(y[j +2]+0.5*k2[j+2]),(t+0.5*d)); else if(j==3) k3[j]=d*f3((y[j-2]+0.5*k2[j-2]),(y[j-1]+0.5*k2[j-1]),(y[j]+0.5*k2[j]),(y[j +1]+0.5*k2[j+1]),(t+0.5*d)); else if(j==4) k3[j]=d*f4((y[j-3]+0.5*k2[j-3]),(y[j-2]+0.5*k2[j-2]),(y[j-1]+0.5*k2[j-1]),(y [j]+0.5*k2[j]),(t+0.5*d)); } for(j=1;j<=m;j++) { if(j==1) k4[j]=d*f1((y[j]+k3[j]),(y[j+1]+k3[j+1]),(y[j+2]+k3[j+2]),(y[j+3]+k3[j+3]),(t +d)); else if(j==2) k4[j]=d*f2((y[j-1]+k3[j-1]),(y[j]+k3[j]),(y[j+1]+k3[j+1]),(y[j+2]+k3[j+2]),(t +d)); else if(j==3) k4[j]=d*f3((y[j-2]+k3[j-2]),(y[j-1]+k3[j-1]),(y[j]+k3[j]),(y[j+1]+k3[j+1]),(t +d)); else if(j==4) k4[j]=d*f4((y[j-3]+k3[j-3]),(y[j-2]+k3[j-2]),(y[j-1]+k3[j-1]),(y[j]+k3[j]),(t +d)); } for(j=1;j<=m;j++) y[j]=y[j]+((k1[j]+2.0*k2[j]+2.0*k3[j]+k4[j])/6.0); t=i*d; Db[i]=X-r[i]; Dc[i]=Y-C[i]; //传递t不断更新预测初值 B[i+1]=B[0]+M[0]*t; C[i+1]=C[0]+M[1]*t; G[i+1]=G[0]+M[2]*t; r[i+1]=sqrt(B[i+1]*B[i+1]+G[i+1]*G[i+1]); R[i+1]=sqrt(r[i+1]*r[i+1]+C[i+1]*C[i+1]); if((0<Db[i])&&(Db[i]<=1))//误差选择1m只能保证一个误差,另一个相差很大 { //if((0<Dc[i])&&(Dc[i]<=1))//加上此判断条件时角度为负值,按理论来说不应该是负值 A[i]=atan(Y/X);//此时为弧度 A[i]=A[i]*180/PI; if(i%1==0) printf("t=%.1f,角度A=%.1f,方程2X=%.1f,Y=%.1f,方程1r=%.1f,C=%.1f\n",t,A[i],X,Y,r[i +1],C[i+1]); } else//重新更改方程2的迭代初值, { A[i+1]=asin(A[i]); A[i+1]=A[i+1]+Dc[i]/R[i]; A[i+1]=A[i+1]*180/PI; A[i+1]=cos(A[i+1]*PI/180); y0[1]=V*A[i+1]; A[i+1]=asin(A[i]); A[i+1]=A[i+1]+Dc[i]/R[i]; A[i+1]=A[i+1]*180/PI; A[i+1]=sin(A[i+1]*PI/180); y0[2]=V*A[i+1]; y0[3]=0; y0[4]=0; for(j=1;j<=m;j++); y[j]=y0[j]; // 加上下面程序(改变步长d)则运行停止工作,没有时运行良好? /* if((0.3*r[i]<X)&&(X<0.6*r[i])) { d=0.7*d; if((0.6*r[i]<X)&&(X<0.9*r[i])) d=0.4*d; if((0.9*r[i]<X)&&(X<r[i])) d=0.1*d; }*/ } } }

利用四阶龙格库塔法求解一阶微分组的c++程序

#include <iostream> #include <iomanip> using namespace std; void RK4( double (*f)(double t,double x, double y),double (*g)(double t,double x, double y) ,double initial[3], double resu[3],double h) { double f1,f2,f3,f4,g1,g2,g3,g4,t0,x0,y0,x1,y1; t0=initial[0];x0=initial[1];y0=initial[2]; f1=f(t0,x0,y0); f2=f(t0+h/2, x0+h*f1/2,y0+h*g1/2); f3=f(t0+h/2, x0+h*f2/2,y0+h*g2/2); f4=f(t0+h, x0+h*f3,y0+h*g3); g1=g(t0,x0,y0); g2=g(t0+h/2,x0+h*f1/2,y0+h*g1/2); g3=g(t0+h/2,x0+h*f2/2,y0+h*g2/2); g4=g(t0+h, x0+h*f3,y0+h*g3); x1=x0+h*(f1+2*f2+2*f3+f4)/6; y1=y0+h*(g1+2*g2+2*g3+g4)/6; resu[0]=t0+h;resu[1]=x1;resu[2]=y1; } int main() { double f(double t,double x, double y); double g(double t,double x, double y); double initial[3],resu[3]; double a,b,H; double t,step; int i; cout<<"输入所求微分方程组的初值t0,x0,y0:"; cin>>initial[0]>>initial[1]>>initial[2]; cout<<"输入所求微分方程组的微分区间[a,b]:"; cin>>a>>b; cout<<"输入所求微分方程组所分解子区间的个数step:"; cin>>step; //cout<<setiosflags(ios::right)<<setiosflags(ios::fixed)<<setprecision(10); H=(b-a)/step; cout<< initial[0]<<setw(18)<<initial[1]<<setw(18)<<initial[2]<<endl; for(i=0;i<step;i++) { RK4( f,g ,initial, resu,H); cout<<resu[0]<<setw(20)<<resu[1]<<setw(20)<<resu[2]<<endl; initial[0]=resu[0];initial[1]=resu[1];initial[2]=resu[2]; } return(0); } double f(double t,double x, double y) { double dx; dx=x*(300-x/10)*(1-x/300-y*20/300); return(dx); } double g(double t,double x, double y) { double dy; dy=y*(300-y/10)*(1-y/300-x*1/20/300); return(dy); } 这是利用四阶龙格库塔法求解一阶微分组的c++程序,但是当我将t0,x0,y0分别设置为0 ,50 ,50,t范围取[0 1000],求解次数为1000次时为啥求解不出来??? 求教各位了,十分着急!希望懂得人为小弟解惑

求教用c++编写的龙格库塔公式,要适应于微分方程组的,不是单个方程,最好有实例,谢谢各路大神!!

求教用c++编写的龙格库塔公式,要适应于微分方程组的,不是单个方程,最好有实例,谢谢各路大神!!

龙哥库塔4阶matlab求初值问题

![图片说明](https://img-ask.csdn.net/upload/201611/21/1479736732_300802.png)![图片说明](https://img-ask.csdn.net/upload/201611/21/1479736776_916346.png)![图片说明](https://img-ask.csdn.net/upload/201611/21/1479736817_424237.png)不知道哪里出错了,一直报错!如下图所示:![图片说明](https://img-ask.csdn.net/upload/201611/21/1479736881_290288.png)希望有路过的大神助攻一下!这个问题困扰了我一晚上了。。。。

关于C语言实现龙格库塔法的问题

能有高手编一个用四阶龙格库塔法求解初值问题(微分方程)的程序吗,有急用?

C++龙格库塔法步长问题

关于用龙哥库塔法求微分方程的问题: 在程序里每次改变步长h,求得结果都会不一样比如取步长为0.2和0.5,计算结果就不一样,不知道是什么原因,请大神帮帮我,谢谢,我把那个微分方程以及代码附上: > 微分方程:dy/dx=y*y;0<=x<=1; y(0)=1; 代码: double f(double x, double y)定义函数f,用来求微分方程 { return y*y; } /////////////////////////// void CQweDlg::OnButton1() { // TODO: Add your control notification handler code here double xx[100][100],yy[100][100],k[100][100]; double h;步长 int i; char ch1[10]; GetDlgItem(IDC_EDIT1)->GetWindowText(ch1,10); h=atof(ch1); long double T; T=1/h;循环次数 xx[1][0]=0; yy[1][0]=1; for(i=0; i<=T;i++) { k[1][1]=f(xx[1][i],yy[1][i]); k[1][2]=f(xx[1][i]+h/2,yy[1][i]+h*k[1][1]/2); k[1][3]=f(xx[1][i]+h/2,yy[1][i]+h*k[1][2]/2); k[1][4]=f(xx[1][i]+h,yy[1][i]+h*k[1][3]); yy[1][i+1]=yy[1][i]+h*(k[1][1]+2*k[1][2]+2*k[1][3]+k[1][4])/6; xx[1][i+1]=xx[1][i]+h; CString str1,str2; str1.Format("%.20f",yy[1][1000]); str2.Format("%.2f",h); GetDlgItem(IDC_EDIT2)->SetWindowText(str1); } }

求随机微分方程解的matlab仿真程序

dx=f(x)dt+g(x)dw,其中w 为布朗运动,求随机微分方程解的matlab仿真程序

图形化编程实现改进的欧拉格式和龙格库塔格式。这里有个C语言的,想改写成C#。

1)改进欧拉法求解常微分方程的初值问题 #include <stdio.h> float func(float x,float y) { return(y-x); } float euler(float x0,float xn,float y0,int N) { float x,y,yp,yc,h; int i; x=x0; y=y0; h=(xn-x0)/(float)N; for(i=1;i<=N;i++) { yp=y+h*func(x,y); x=x0+i*h; yc=y+h*func(x,yp); y=(yp+yc)/2.0; } return(y); } main() { float x0,xn,y0,e; int n; printf("\ninput n:\n "); scanf("%d",&n); printf("input x0,xn:\n "); scanf("%f,%f",&x0,&xn); printf("input y0:\n "); scanf("%f",&y0); e=euler(x0,xn,y0,n); printf("y(%f)=%6.4f",y0,e); } input n: 20 input x0,xn: 1,6 input y0: 2 y(2.000000)=7.0000Press any key to continue (2)四阶龙格—库塔法 #include <stdio.h> float func(float x,float y) { return(x-y); } float runge_kutta(float x0,float xn,float y0,int N) { float x,y,y1,y2,h,xh; float d1,d2,d3,d4; int i; x=x0; y=y0; h=(xn-x0)/(float)N; for(i=1;i<=N;i++) { xh=x+h/2; d1=func(x,y); d2=func(xh,y+h*d1/2.0); d3=func(xh,y+h*d2/2.0); d4=func(xh,y+h*d3); y=y+h*(d1+2*d2+2*d3+d4)/6.0; x=x0+i*h; } return(y); } main() { float x0,xn,y0,e; int N; printf("\ninput n:\n "); scanf("%d",&N); printf("input x0,xn:\n "); scanf("%f,%f",&x0,&xn); printf("input y0:\n "); scanf("%f",&y0); e=runge_kutta(x0,xn,y0,N); printf("y(%f)=%8.6f",y0,e); } input n: 10 input x0,xn: 1,2 input y0: 5 y(5.000000)=2.833863Press any key to continue

如何求解Lotka-Volterra竞争方程

求教如何c++编程求解Lotka-Volterra竞争方程?

MATLAB求解多元高阶微分方程组

x,y,z为位移,一阶导为速度,二阶导为加速度,m,F,c和两个角度为常数,用MATLAB求解微分方程组,![求程序e](e\1.png)

runge-kutta法求解常微分方程组

x(t)'=(t-1)/0.001 I2(t)'=cos(x) 初值为:x(0)=0 I2(0)=1 这种套的常微分方程如何用runge-kutta法求解,直接套用runge-kutta法求解出来的图像是错的。各位程序界大神们帮帮忙。谢谢大家

求一个VC下实现的C++代码

曲线绘制 根据以下微分方程可以产生任意多的实时数据: 其中参数 ,初始条件取 ,且当 时。用4阶龙格库塔方法对上述方程实现离散化。 要求显示实时数据曲线,并能用左右箭头、PgUp、PgDn实现曲线的滚动和翻页。

MATLAB用euler方法和方向场刻画微分方程的解,如何理解呢

微分方程为:dx/dt=m.*A(P-(n.*R.*T)/(V-x)) 其中m=0.1、A=1.5、P=15、(n.*R.*T)=9.6814、V=4 如何用euler方法求解这个微分方程呢,如何用方向场来刻画呢,使用MATLAB来实现 求代码,感谢大佬

matlab赋值时元素个数不同

function [minn, maxx] = bianjie(h1) %guiyihua中调用,从边界图像中,内切提取内边界的值 m = []; [r, c] = find(bwlabel(h1) == 1); m(1) = max(r); m(2) = min(r); [r, c] = find(bwlabel(h1) == 2); m(3) = max(r); m(4) = min(r); if m(1) < m(3) minn = m(1); maxx = m(4); else minn = m(3); maxx = m(2); end 错误原因为:在赋值A(:)=B中,A和B的元素数目必须相同

求助:循环角度不断改变初值,解常微分方程

四阶龙格库塔解微分方程, #include<stdio.h> #include <math.h> #include<stdlib.h> #define PI 3.141593 #define Vx y[1] #define Vy y[2] #define X y[3] #define Y y[4] //方程组2 #define f1(y1,y2,y3,y4,t) (-0.1109*y1)//Vx(Vr) #define f2(y1,y2,y3,y4,t) (-0.1109*y2-9.81)//Vy(Vh) #define f3(y1,y2,y3,y4,t) (y1)//X(r) #define f4(y1,y2,y3,y4,t) (y2)//Y(h) void main(void) { int k,s,i,j,V,m=4,n=150; double Db[151],Dc[151],T,t,a[151],A[151],A1[151],B[151],C[151],G[151],M [151],r,R,d,y[151],y0[151],k1[151],k2[151],k3[151],k4[151]; B[0]=0; C[0]=100; G[0]=100; M[0]=0; M[1]=0; M[2]=-10; V=120; d=0.1; T=0; for(s=1;s<=100;s++) { B[s]=B[s-1]+M[0]*T;//T时刻方程组1 C[s]=C[s-1]+M[1]*T; G[s]=G[s-1]+M[2]*T; r=sqrt(B[s]*B[s]+G[s]*G[s]); R=sqrt(r*r+C[s]*C[s]); A[s]=atan(C[s]/r); A[s]=A[s]*180/PI; A[s]=cos(A[s]*PI/180);//输入必须为弧度 y0[1]=V*A[s]; A[s]=atan(C[s]/r); A[s]=A[s]*180/PI;//弧度化为度 A[s]=sin(A[s]*PI/180);//输入必须为弧度 y0[2]=V*A[s]; y0[3]=0; y0[4]=0; t=0; for(j=1;j<=m;j++) y[j]=y0[j]; for(i=1;i<n;i++) { for(j=1;j<=m;j++) { if(j==1) k1[j]=d*f1((y[j]),(y[j+1]),(y[j+2]),(y[j+3]),t); else if(j==2) k1[j]=d*f2((y[j-1]),(y[j]),(y[j+1]),(y[j+2]),t); else if(j==3) k1[j]=d*f3((y[j-2]),(y[j-1]),(y[j]),(y[j+1]),t); else if(j==4) k1[j]=d*f4((y[j-3]),(y[j-2]),(y[j-1]),(y[j]),t); } for(j=1;j<=m;j++) { if(j==1) k2[j]=d*f1((y[j]+0.5*k1[j]),(y[j+1]+0.5*k1[j+1]),(y[j+2]+0.5*k1[j+2]),(y[j +3]+0.5*k1[j+3]),(t+0.5*d)); else if(j==2) k2[j]=d*f2((y[j-1]+0.5*k1[j-1]),(y[j]+0.5*k1[j]),(y[j+1]+0.5*k1[j+1]),(y[j +2]+0.5*k1[j+2]),(t+0.5*d)); else if(j==3) k2[j]=d*f3((y[j-2]+0.5*k1[j-2]),(y[j-1]+0.5*k1[j-1]),(y[j]+0.5*k1[j]),(y[j +1]+0.5*k1[j+1]),(t+0.5*d)); else if(j==4) k2[j]=d*f4((y[j-3]+0.5*k1[j-3]),(y[j-2]+0.5*k1[j-2]),(y[j-1]+0.5*k1[j-1]), (y[j]+0.5*k1[j]),(t+0.5*d)); } for(j=1;j<=m;j++) { if(j==1) k3[j]=d*f1((y[j]+0.5*k2[j]),(y[j+1]+0.5*k2[j+1]),(y[j+2]+0.5*k2[j+2]),(y[j +3]+0.5*k2[j+3]),(t+0.5*d)); else if(j==2) k3[j]=d*f2((y[j-1]+0.5*k2[j-1]),(y[j]+0.5*k2[j]),(y[j+1]+0.5*k2[j+1]),(y[j +2]+0.5*k2[j+2]),(t+0.5*d)); else if(j==3) k3[j]=d*f3((y[j-2]+0.5*k2[j-2]),(y[j-1]+0.5*k2[j-1]),(y[j]+0.5*k2[j]),(y[j +1]+0.5*k2[j+1]),(t+0.5*d)); else if(j==4) k3[j]=d*f4((y[j-3]+0.5*k2[j-3]),(y[j-2]+0.5*k2[j-2]),(y[j-1]+0.5*k2[j-1]), (y[j]+0.5*k2[j]),(t+0.5*d)); } for(j=1;j<=m;j++) { if(j==1) k4[j]=d*f1((y[j]+k3[j]),(y[j+1]+k3[j+1]),(y[j+2]+k3[j+2]),(y[j+3]+k3[j+3]), (t+d)); else if(j==2) k4[j]=d*f2((y[j-1]+k3[j-1]),(y[j]+k3[j]),(y[j+1]+k3[j+1]),(y[j+2]+k3[j+2]), (t+d)); else if(j==3) k4[j]=d*f3((y[j-2]+k3[j-2]),(y[j-1]+k3[j-1]),(y[j]+k3[j]),(y[j+1]+k3[j+1]), (t+d)); else if(j==4) k4[j]=d*f4((y[j-3]+k3[j-3]),(y[j-2]+k3[j-2]),(y[j-1]+k3[j-1]),(y[j]+k3[j]), (t+d)); } for(j=1;j<=m;j++) y[j]=y[j]+((k1[j]+2.0*k2[j]+2.0*k3[j]+k4[j])/6.0); t=i*d; Db[i]=X-r; Dc[i]=Y-C[s];//误差选择5m,方程组2的值X、Y需要不断接近方程组1解值C、G,达到精度输出此 时t、A if(0<Dc[i]<=5) { A[i]=asin(A[i]);//此时为弧度 A[i]=A[i]*180/PI; if(i%10==0) printf("时间=%.2lf,射角A=%.1lf,击中点X=%.1lf,Y=%.1lf\n",t,A[i],X,Y); } else//修正A,更新迭代初值 { A[s]=atan(C[s]/r); A[s]=A[s]+Db[i]/R; A[s]=A[s]*180/PI; A[s]=cos(A[s]*PI/180); y0[1]=V*A[s]; A[s]=atan(C[s]/r);//弹丸射角 A[s]=A[s]+Db[i]/R; A[s]=A[s]*180/PI; A[s]=sin(A[s]*PI/180); y0[2]=V*A[s]; y0[3]=0; y0[4]=0; for(j=1;j<=m;j++); y[j]=y0[j]; } } //超过最大迭代次数时,重新解方程组1 T=T+n*d; } }

关于九轴陀螺仪的数据融合问题

最近在弄九轴陀螺仪的数据融合,用的AHRS算法,获取数据后,如果不加磁力计数据,roll轴和pitch轴数据都还正常,yaw轴会一直偏移,可是加上磁力计后数据就不正常了,会响应很慢,而且还各个轴输出数据都不准确。用的就是网上开源的数据融合代码。想问问可能是哪些方面的问题? ``` /* Feed the sensor data in NED(North East Down) reference frame. Rotation can be added. */ int AHRS_Update2(fp32 QUAT[4], fp32 ACC[3], fp32 GYRO[3], fp32 MAG[3]) { float recipNorm; float q0q0, q0q1, q0q2, q0q3, q1q1, q1q2, q1q3, q2q2, q2q3, q3q3; float hx, hy, hz, bx, bz; float halfvx, halfvy, halfvz, halfwx, halfwy, halfwz; float halfex, halfey, halfez; float qa, qb, qc; if (MAG == NULL) return AHRS_UpdateIMU(QUAT, ACC, GYRO); float mx = MAG[0]; float my = MAG[1]; float mz = MAG[2]; // Use IMU algorithm if magnetometer measurement invalid (avoids NaN in magnetometer normalisation) if((mx == 0.0f) && (my == 0.0f) && (mz == 0.0f)) { return AHRS_UpdateIMU(QUAT, ACC, GYRO); } float ax = ACC[0]; float ay = ACC[1]; float az = ACC[2]; float gx = GYRO[0]; float gy = GYRO[1]; float gz = GYRO[2]; // Compute feedback only if accelerometer measurement valid (avoids NaN in accelerometer normalisation) // 只在加速度计数据有效时才进行运算 if(!((ax == 0.0f) && (ay == 0.0f) && (az == 0.0f))) { // Normalise accelerometer measurement // 将加速度计得到的实际重力加速度向量v单位化 recipNorm = invSqrt(ax * ax + ay * ay + az * az); ax *= recipNorm; ay *= recipNorm; az *= recipNorm; // Normalise magnetometer measurement // 将磁力计得到的实际磁场向量m单位化 recipNorm = invSqrt(mx * mx + my * my + mz * mz); mx *= recipNorm; my *= recipNorm; mz *= recipNorm; // Auxiliary variables to avoid repeated arithmetic // 辅助变量,以避免重复运算 q0q0 = QUAT[0] * QUAT[0]; q0q1 = QUAT[0] * QUAT[1]; q0q2 = QUAT[0] * QUAT[2]; q0q3 = QUAT[0] * QUAT[3]; q1q1 = QUAT[1] * QUAT[1]; q1q2 = QUAT[1] * QUAT[2]; q1q3 = QUAT[1] * QUAT[3]; q2q2 = QUAT[2] * QUAT[2]; q2q3 = QUAT[2] * QUAT[3]; q3q3 = QUAT[3] * QUAT[3]; // Reference direction of Earth's magnetic field // 通过磁力计数据与坐标转换矩阵得到理论磁场向量 // 公式(5) hx = 2.0f * (mx * (0.5f - q2q2 - q3q3) + my * (q1q2 - q0q3) + mz * (q1q3 + q0q2)); hy = 2.0f * (mx * (q1q2 + q0q3) + my * (0.5f - q1q1 - q3q3) + mz * (q2q3 - q0q1)); hz = 2.0f * (mx * (q1q3 - q0q2) + my * (q2q3 + q0q1) + mz * (0.5f - q1q1 - q2q2)); // 公式(6) bx = sqrt(hx * hx + hy * hy); bz = 2.0f * (mx * (q1q3 - q0q2) + my * (q2q3 + q0q1) + mz * (0.5f - q1q1 - q2q2)); // Estimated direction of gravity and magnetic field // 将理论重力加速度向量与理论磁场向量通过坐标转换矩阵转换至机体坐标系 // 注意,这里实际上是矩阵*1/2,在开头对Kp Ki的宏定义均为2*增益 // 至于这么设计程序的原因笔者也没有弄清,猜测可能是为减少乘法运算量? // 笔记(一)中内容 halfvx = q1q3 - q0q2; halfvy = q0q1 + q2q3; halfvz = q0q0 - 0.5f + q3q3; // 公式(7) halfwx = bx * (0.5f - q2q2 - q3q3) + bz * (q1q3 - q0q2); halfwy = bx * (q1q2 - q0q3) + bz * (q0q1 + q2q3); halfwz = bx * (q0q2 + q1q3) + bz * (0.5f - q1q1 - q2q2); // Error is sum of cross product between estimated direction and measured direction of field vectors // 通过向量外积得到重力加速度向量和磁场向量的实际值与测量值直接误差 // 公式(8)(9)(10) halfex = (ay * halfvz - az * halfvy) + (my * halfwz - mz * halfwy); halfey = (az * halfvx - ax * halfvz) + (mz * halfwx - mx * halfwz); halfez = (ax * halfvy - ay * halfvx) + (mx * halfwy - my * halfwx); // Compute and apply integral feedback if enabled // 在误差补偿PI控制器中积分项使能情况下计算并应用积分项 if(twoKi > 0.0f) { // integral error scaled by Ki // 积分过程 integralFBx += twoKi * halfex * (1.0f / sampleFreq); integralFBy += twoKi * halfey * (1.0f / sampleFreq); integralFBz += twoKi * halfez * (1.0f / sampleFreq); // apply integral feedback // 应用误差补偿中的积分项 gx += integralFBx; gy += integralFBy; gz += integralFBz; } else { // prevent integral windup // 避免为负值的Ki时积分异常饱和 integralFBx = 0.0f; integralFBy = 0.0f; integralFBz = 0.0f; } // Apply proportional feedback // 应用误差补偿中的比例项 gx += twoKp * halfex; gy += twoKp * halfey; gz += twoKp * halfez; } // Integrate rate of change of quaternion //一阶龙格库塔法迭代四元数 gx *= (0.5f * (1.0f / sampleFreq)); // pre-multiply common factors gy *= (0.5f * (1.0f / sampleFreq)); gz *= (0.5f * (1.0f / sampleFreq)); qa = QUAT[0]; qb = QUAT[1]; qc = QUAT[2]; QUAT[0] += (-qb * gx - qc * gy - QUAT[3] * gz); QUAT[1] += (qa * gx + qc * gz - QUAT[3] * gy); QUAT[2] += (qa * gy - qb * gz + QUAT[3] * gx); QUAT[3] += (qa * gz + qb * gy - qc * gx); // Normalise quaternion // 单位化四元数 保证四元数在迭代过程中保持单位性质 recipNorm = invSqrt(QUAT[0] * QUAT[0] + QUAT[1] * QUAT[1] + QUAT[2] * QUAT[2] + QUAT[3] * QUAT[3]); QUAT[0] *= recipNorm; QUAT[1] *= recipNorm; QUAT[2] *= recipNorm; QUAT[3] *= recipNorm; return 0; } ```

关于mpu9250原始数据的滤波

因为加速度计和磁力计存在高频误差,需要低通滤波器来输出准确的数据,陀螺仪存在低频误差,需要进行高通滤波,但是我看了CSDN上的帖子,也没帖子明确这些滤波的截止频率,所以想来这请教一下各位。 ![图片说明](https://img-ask.csdn.net/upload/201910/28/1572232194_532982.png) 这个是我在帖子上看到的关于mpu9250解算姿态的流程,我最后的姿态融合用的是AHRS算法,进行了磁力计的补偿,但是不知道这种情况下还需要进行第二步(磁力计坐标系到加速度坐标系)吗?需要的画,相关的资料哪位大哥有啊,有的话能提供一下吗?救救老弟吧,真的谢谢了。 下面是我用的AHRS算法代码: ``` // 加速度计、地磁计、陀螺仪数据融合,更新四元数 /* [gx,gy,gz]为陀螺仪的测量值 [ax,at,az]为加速度的测量值 [mx,my,mz]为地磁计的测量值 */ void AHRSupdate(float gx, float gy, float gz, float ax, float ay, float az, float mx, float my, float mz) { float norm; float hx, hy, hz, bx, bz; float vx, vy, vz, wx, wy, wz; float ex, ey, ez; // 定义一些辅助变量用于转换矩阵 float q0q0 = q0*q0; float q0q1 = q0*q1; float q0q2 = q0*q2; float q0q3 = q0*q3; float q1q1 = q1*q1; float q1q2 = q1*q2; float q1q3 = q1*q3; float q2q2 = q2*q2; float q2q3 = q2*q3; float q3q3 = q3*q3; // 归一化加速度计和地磁计的度数 norm = sqrt(ax*ax + ay*ay + az*az); ax = ax / norm; ay = ay / norm; az = az / norm; norm = sqrt(mx*mx + my*my + mz*mz); mx = mx / norm; my = my / norm; mz = mz / norm; //将b系中的地磁计分量[mx,my,mz]转换到n系,得到[hx,hy,hz] hx = 2*mx*(0.5 - q2q2 - q3q3) + 2*my*(q1q2 - q0q3) + 2*mz*(q1q3 + q0q2); hy = 2*mx*(q1q2 + q0q3) + 2*my*(0.5 - q1q1 - q3q3) + 2*mz*(q2q3 - q0q1); hz = 2*mx*(q1q3 - q0q2) + 2*my*(q2q3 + q0q1) + 2*mz*(0.5 - q1q1 - q2q2); //得到n系中的地磁向量的真实值[bx,bz,by],其中by=0 bx = sqrt((hx*hx) + (hy*hy)); bz = hz; //n系中的地磁向量[bx,by,bz]转换到b系中,得到[wx,wy,wz] wx = 2*bx*(0.5 - q2q2 - q3q3) + 2*bz*(q1q3 - q0q2); wy = 2*bx*(q1q2 - q0q3) + 2*bz*(q0q1 + q2q3); wz = 2*bx*(q0q2 + q1q3) + 2*bz*(0.5 - q1q1 - q2q2); //n系中重力加速度[0,0,1]转换到b系中得到三个分量[vx,vy,vz] vx = 2*(q1q3 - q0q2); vy = 2*(q0q1 + q2q3); vz = q0q0 - q1q1 - q2q2 + q3q3; //计算[wx,wy,wz] X [mx,my,mz],[ax,at,az] X [vx,vy,vz],得到两个误差后求和 ex = (ay*vz - az*vy) + (my*wz - mz*wy); ey = (az*vx - ax*vz) + (mz*wx - mx*wz); ez = (ax*vy - ay*vx) + (mx*wy - my*wx); //PI控制器中的积分部分 exInt = exInt + ex*Ki* (1.0f / sampleFreq); eyInt = eyInt + ey*Ki* (1.0f / sampleFreq); ezInt = ezInt + ez*Ki* (1.0f / sampleFreq); //误差经过PI控制器后输出,然后补偿到角速度的三个分量,Kp、Ki是需要调节的参数 gx = gx + Kp*ex + exInt; gy = gy + Kp*ey + eyInt; gz = gz + Kp*ez + ezInt; //一阶龙格库塔法更新四元数 q0 = q0 + (-q1*gx - q2*gy - q3*gz)*halfT; q1 = q1 + (q0*gx + q2*gz - q3*gy)*halfT; q2 = q2 + (q0*gy - q1*gz + q3*gx)*halfT; q3 = q3 + (q0*gz + q1*gy - q2*gx)*halfT; // 归一化四元数 norm = sqrt(q0*q0 + q1*q1 + q2*q2 + q3*q3); q0 = q0 / norm; q1 = q1 / norm; q2 = q2 / norm; q3 = q3 / norm; } ```

YOLOv3目标检测实战:训练自己的数据集

YOLOv3目标检测实战:训练自己的数据集

150讲轻松搞定Python网络爬虫

150讲轻松搞定Python网络爬虫

实用主义学Python(小白也容易上手的Python实用案例)

实用主义学Python(小白也容易上手的Python实用案例)

我说我不会算法,阿里把我挂了。

不说了,字节跳动也反手把我挂了。

立方体线框模型透视投影 (计算机图形学实验)

计算机图形学实验 立方体线框模型透视投影 的可执行文件,亲测可运行,若需报告可以联系我,期待和各位交流

2019 AI开发者大会

2019 AI开发者大会

组成原理课程设计(实现机器数的真值还原等功能)

实现机器数的真值还原(定点小数)、定点小数的单符号位补码加减运算、定点小数的补码乘法运算和浮点数的加减运算。

C/C++跨平台研发从基础到高阶实战系列套餐

一 专题从基础的C语言核心到c++ 和stl完成基础强化; 二 再到数据结构,设计模式完成专业计算机技能强化; 三 通过跨平台网络编程,linux编程,qt界面编程,mfc编程,windows编程,c++与lua联合编程来完成应用强化 四 最后通过基于ffmpeg的音视频播放器,直播推流,屏幕录像,

MFC一站式终极全套课程包

该套餐共包含从C小白到C++到MFC的全部课程,整套学下来绝对成为一名C++大牛!!!

软件测试2小时入门

软件测试2小时入门

三个项目玩转深度学习(附1G源码)

三个项目玩转深度学习(附1G源码)

计算机图形学-球的光照模型课程设计

计算机图形学-球的光照模型,有代码完美运行,有课程设计书

Linux常用命令大全(非常全!!!)

Linux常用命令大全(非常全!!!) 最近都在和Linux打交道,感觉还不错。我觉得Linux相比windows比较麻烦的就是很多东西都要用命令来控制,当然,这也是很多人喜欢linux的原因,比较短小但却功能强大。我将我了解到的命令列举一下,仅供大家参考: 系统信息 arch 显示机器的处理器架构 uname -m 显示机器的处理器架构 uname -r 显示正在使用的内核版本 d...

因为看了这些书,我大二就拿了华为Offer

四年了,四年,你知道大学这四年我怎么过的么?

深度学习原理+项目实战+算法详解+主流框架(套餐)

深度学习系列课程从深度学习基础知识点开始讲解一步步进入神经网络的世界再到卷积和递归神经网络,详解各大经典网络架构。实战部分选择当下最火爆深度学习框架PyTorch与Tensorflow/Keras,全程实战演示框架核心使用与建模方法。项目实战部分选择计算机视觉与自然语言处理领域经典项目,从零开始详解算法原理,debug模式逐行代码解读。适合准备就业和转行的同学们加入学习! 建议按照下列课程顺序来进行学习 (1)掌握深度学习必备经典网络架构 (2)深度框架实战方法 (3)计算机视觉与自然语言处理项目实战。(按照课程排列顺序即可)

fakeLocation13.5.1.zip

fakeLocation13.5.1 虚拟定位 ios13.5.1的最新驱动下载,iPhone/iPad免越狱虚拟定位工具Location-cleaned驱动已更新

UnityLicence

UnityLicence

Python可以这样学(第一季:Python内功修炼)

Python可以这样学(第一季:Python内功修炼)

Python+OpenCV计算机视觉

Python+OpenCV计算机视觉

土豆浏览器

土豆浏览器可以用来看各种搞笑、电影、电视剧视频

【数据结构与算法综合实验】欢乐连连看(C++ & MFC)案例

这是武汉理工大学计算机学院数据结构与算法综合实验课程的第三次项目:欢乐连连看(C++ & MFC)迭代开发代码。运行环境:VS2017。已经实现功能:开始游戏、消子、判断胜负、提示、重排、计时、帮助。

php+mysql学生成绩管理系统

学生成绩管理系统,分三个模块:学生,教师和管理员。 管理员模块:负责学生、老师信息的增删改;发布课程信息的增删改,以便让学生选课;审核老师提交的学生成绩并且打印成绩存档;按照课号查询每个课号的学生成绩

多功能数字钟.zip

利用数字电子计数知识设计并制作的数字电子钟(含multisim仿真),该数字钟具有显示星期、24小时制时间、闹铃、整点报时、时间校准功能

推荐24个国外黄色网站欣赏

在中国清朝,明黄色的衣服只有皇子才有资格穿,慢慢的黄色在中国就成了高贵的颜色。在人们的色彩印象中,黄色也表现为暂停。所以当你的网页设计采用黄色的时候,会让人们在你的网页前停留。 黄色,就像橙色和红色,黄色也是一个暖色。它有大自然、阳光、春天的涵义,而且通常被认为是一个快乐和有希望的色彩。黄色是所有色相中最能发光的颜色,给人轻快,透明,辉煌,充满希望的色彩印象。 黄色是一个高可见的色...

u-boot-2015.07.tar.bz2

uboot-2015-07最新代码,喜欢的朋友请拿去

一学即懂的计算机视觉(第一季)

一学即懂的计算机视觉(第一季)

学生成绩管理系统(PHP + MYSQL)

做的是数据库课程设计,使用的php + MySQL,本来是黄金搭配也就没啥说的,推荐使用wamp服务器,里面有详细的使用说明,带有界面的啊!呵呵 不行的话,可以给我留言!

Windows版YOLOv4目标检测实战:训练自己的数据集

Windows版YOLOv4目标检测实战:训练自己的数据集

C++语言基础视频教程

C++语言基础视频教程

玩转Python-Python3基础入门

玩转Python-Python3基础入门

2019校招硬件乐鑫+比特大陆笔试题

楼主水硕一枚,参加了2019年的秋招。自己总结了下乐鑫的笔试题目(现场笔试)以及网上考试的比特大陆的题目

相关热词 c# 开发接口 c# 中方法上面的限制 c# java 时间戳 c#单元测试入门 c# 数组转化成文本 c#实体类主外键关系设置 c# 子函数 局部 c#窗口位置设置 c# list 查询 c# 事件 执行顺序
立即提问
相关内容推荐