troubles with ode45 - Programmers Heaven

Howdy, Stranger!

It looks like you're new here. If you want to get involved, click one of these buttons!

Categories

Welcome to the new platform of Programmer's Heaven! We apologize for the inconvenience caused, if you visited us from a broken link of the previous version. The main reason to move to a new platform is to provide more effective and collaborative experience to you all. Please feel free to experience the new platform and use its exciting features. Contact us for any issue that you need to get clarified. We are more than happy to help you.

troubles with ode45

NotYetMScNotYetMSc Posts: 1Member
Hello everybody!
I've written a program:
[code]
clear all
close all
clc

global mk Kfi g lk alfa beta Kz Jk a0 nsteps1
g=9.81;
E=7*10^10;
alfa=20*pi/180;
L=50;
N=5;
mk=300;
lk=1.5;
H=40000;
ro1=8.45;
d=0.04;
v0=0;
a0=0.5;
%**************************************************************************
roz1=ro1+(N*mk)/L;
Jk=mk*lk^2;
I=d^4*pi/64;
Kfi=(L/(N*cos(alfa)))/6*E*I;
%**************************************************************************
h=0.01; %step
%**************************************************************************
droga=L/cos(alfa);
delta=v0^2+4*(a0/2)*droga;
tc=(-v0+sqrt(delta))/a0;
Vmax=3;
V=v0+a0*tc;
if V>Vmax
t1=(Vmax-v0)/a0;
s1=v0*t1+((a0*t1^2)/2);
s2=droga-s1;
t2=s2/Vmax;
tc=t1+t2;
t01=[0:h:t1];
t02=[h:h:t2];
Lk1=cos(alfa).*((v0.*t01+(a0.*t01.^2))./2);
Lk2=Lk1(1,end)+(cos(alfa).*(Vmax.*t02));
Lk=[Lk1 Lk2];
t0=[0:h:tc];
%modyfying (x/0)
Lkk=Lk;
Lkk(Lkk<1)=1;
Lkk(Lkk>(L-1))=(L-1);
z=((roz1*g)/(2*H*(cos(alfa))^2)).*Lk.*(L-Lk);
beta=atan((roz1*g/H)*(L/2-Lk));
prosta=-tan(alfa).*Lk;
tor=prosta-z;
Kz=H*((1./Lkk)+(1./(L-Lkk)));
%beta=beta';
%beta=median(beta,1);
%Kz=Kz';
%Kz=median(Kz,1);
nsteps1=t1/h;
%nsteps2=t2/h;
tspan1 = linspace(0,t1,nsteps1);
tspan2 = linspace(t1,t2,nsteps2);
[tc1,x1]=ode45(@przysp1, tspan1, [0.02 0 0.0028 0]); [code/]
and function:
[code]
function xprime=przysp1(tc1,x1)
global mk Kfi g lk alfa beta Kz Jk a0 nsteps1
i=1;

xprime=zeros(4,1);
xprime(1)=x1(2);
xprime(2)=(mk*(Kfi*x1(1)+g*mk*lk*sin(x1(1)))*cos(2*(alfa+beta(:,i)))+(x1(2)^2)*(mk^2)*...
(lk^2)*(sin(2*(x1(1)))*cos(2*(alfa+beta(:,i))))-cos(x1(1)-alfa-beta(:,i))*...
sin(x1(1)+alfa+beta(:,i))*a0*(mk^2)*lk*sin(-beta(:,i))*sin(x1(1)+alfa+beta(:,i))-...
cos(alfa-x1(1))*cos(2*(alfa+beta(:,i)))-Kz(:,i)*x1(3)*mk*lk*sin(x1(1)+alfa+beta(:,i)))/...
(mk*cos(2*(alfa+beta(:,i)))*(mk*(lk^2)+Jk)-(mk^2)*(lk^2)*sin(x1(1)+alfa+beta(:,i))*...
sin(x1(1)-alfa-beta(:,i)));
xprime(3)=x1(4);
xprime(4)=(xprime(2)*mk*lk*sin(x1(1)-alfa-beta(:,i))-Kz(:,i)*x1(3)+mk*a0*sin(-beta(:,i))-mk*lk*...
(x1(2)^2)*cos(x1(1)-alfa-beta(:,i)))/(mk*cos(2*(alfa+beta(:,i))));
[code/]

I want that "i" will increase in each step. (in first step i=1, in secound i=2... and so on).

Now program works only with all this values. If I change a0 lk, mk... matlab show me: "Warning: Failure at t=2.802832e+000. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (9.957659e-015) at time t.
(Type "warning off MATLAB:ode45:IntegrationTolNotMet" to suppress this warning.)
> In C:Matlab6p5 oolboxmatlabfunfunode45.m at line 335
In C:Documents and Settings...proba.m at line 64 "
I dont know what is wrong. Please help me... and sorry for my english
Sign In or Register to comment.