error updating function line

i have two problems with this code. first when i run it it solves theta for every step and show me even i don't need to see the exact results of each step. second it gives me error updating function line. how can i fix this?

clear;
clc;
k=3500;
m=5;
c=270;
a=0.75;
T=1.5;
w=2*pi/T;
g=9.81;

syms t

a1=2734;
a2=-6718;
a3=5682;
a4=-1931;
a5=228;
Ft1=100*(a1*t^5+a2*t^4+a3*t^3+a4*t^2+a5*t);
Ft2=-500;

syms n O(i)
a0=1/T * (int(Ft1, t, 0, 1)+int(Ft2, t, 1, 1.5));
an=2/T * (int(Ft1*cos(w*n*t),t,0,1)+ int(Ft2*cos(w*n*t),t,1,1.5));
bn=2/T * (int(Ft1*sin(w*n*t),t,0,1)+ int(Ft2*sin(w*n*t),t,1,1.5));


for i=1:10
    
    
    F(i)=an*cos(w*n*t)+bn*sin(w*n*t);
    Fsum=symsum(F(i),n,1,20);
    Ftoplam= a0 + sum(Fsum);
     theta=O(i)
     4/3*a^2*m* diff(diff(O(i)))+c*a^2*diff(O(i))+k*a^2*O(i)+m*g*a==2*Ftoplam

    if i==10
        break
    end
  i=i+1;
end
 

fplot(O(i),[1 10]);

[Dyuman Joshi]

As you are using a for loop, you don't need to use break or update "i" in each iteration. Remove those lines.

What are these lines supposed to do -

theta = O(i)
4/3*a^2*m* diff(diff(O(i)))+c*a^2*diff(O(i))+k*a^2*O(i)+m*g*a==2*Ftoplam

Do you want to solve the the equation for O and plot it? If yes, then what is the pupose of using theta?

[ege]

there is a theta angle in question which is time dependant. i define it with O(i) to make it function. bottom line is a vibration equation. I need to find O(i) function with 'Ftoplam' values that i calculated previous part.For example lets say 0.4 th second, first i need to find 'Ftoplam' value. Then with that value, i need to calculate O(i) function value at the 0.4th second. after that i need to do it with 0.5th-0.6th..... respectively. final part i need to get time O(i) graph. Hopefully i explained right.

[Torsten]

Do you mean something like this ?

But you will get problems to determine solutions for your 2nd order differential equations for O analytically. Better use a numerical solver by giving two additional boundary conditions for O.

clear;
clc;
k=3500;
m=5;
c=270;
a=0.75;
T=1.5;
w=2*pi/T;
g=9.81;
syms t O(t)
a1=2734;
a2=-6718;
a3=5682;
a4=-1931;
a5=228;
Ft1=100*(a1*t^5+a2*t^4+a3*t^3+a4*t^2+a5*t);
Ft2=-500;
syms n 
a0=1/T * (int(Ft1, t, 0, 1)+int(Ft2, t, 1, 1.5));
a0 = simplify(a0);
an=2/T * (int(Ft1*cos(w*n*t),t,0,1)+ int(Ft2*cos(w*n*t),t,1,1.5));
an = simplify(an);
bn=2/T * (int(Ft1*sin(w*n*t),t,0,1)+ int(Ft2*sin(w*n*t),t,1,1.5));
bn = simplify(bn);

N = 3;
F = cell(1,N);
%sol = cell(1,N);
for i = 1:N
    an_subs = subs(an,n,1:i);
    bn_subs = subs(bn,n,1:i);
    F{i} = a0+sum(an_subs.*cos(w*(1:i)*t)+bn_subs.*sin(w*(1:i)*t));
    %sol{i} = dsolve(4/3*a^2*m* diff(diff(O))+c*a^2*diff(O)+k*a^2*O+m*g*a==2*F{i});
end
fplot(F,[0 10])

%sol

[ege]

O(0) =0

O'(0)=0 and when i wrote this code F results diffeerently than this. Ftoplam should results like in this code.

clear;
clc;
k=3500;
m=5;
c=270;
a=0.75;
T=1.5;
w=2*pi/T;

syms t

a1=2734;
a2=-6718;
a3=5682;
a4=-1931;
a5=228;
Ft1=100*(a1*t^5+a2*t^4+a3*t^3+a4*t^2+a5*t);
Ft2=-500;

syms n
a0=1/T * (int(Ft1, t, 0, 1)+int(Ft2, t, 1, 1.5));
an=2/T * (int(Ft1*cos(w*n*t),t,0,1)+ int(Ft2*cos(w*n*t),t,1,1.5));
bn=2/T * (int(Ft1*sin(w*n*t),t,0,1)+ int(Ft2*sin(w*n*t),t,1,1.5));


for i=1:30
    
     
    F(i)=an*cos(w*n*t)+bn*sin(w*n*t);
    Fsum=symsum(F(i),n,1,20);
    Ftoplam= a0 + sum(Fsum);
    
    if i==30
        break
    end
  i=i+1;
end


fplot(Ftoplam,[0 10])

[Torsten]

What is your loop over i from above good for ? It will always produce the same F(i), the same Fsum and the same Ftoplam.

As far as I understood your problem, you want to solve the 2nd order differential equation

4/3*a^2*m* diff(diff(O))+c*a^2*diff(O)+k*a^2*O+m*g*a==2*F_n

for O for F_n being

F_n(t) = a0 + sum_{j=1}^{j=n} (aj*cos(w*j*t) + bj*sin(w*j*t))

That's what my modification to your code does. Well, it doesn't solve the differential equation because "dsolve" is not able to. But it produces the F_n of the right-hand side and plots them (see above).

If you want something different, please explain it mathematically, not using physical explanations or MATLAB code.

And don't use the i = i+1 command in the for-loop. A for-loop increments i automatically.

Sorry, I forgot the sum in the computation of F:

F = a0+sum(an_subs.*cos(w*(1:i)*t)+bn_subs.*sin(w*(1:i)*t));

instead of

F = a0+an_subs.*cos(w*(1:i)*t)+bn_subs.*sin(w*(1:i)*t);

I modified it in the code above.

Answer 1

[Torsten]

This code plots O(i):

k=3500;
m=5;
c=270;
a=0.75;
T=1.5;
w=2*pi/T;
g=9.81;
syms t 
a1=2734;
a2=-6718;
a3=5682;
a4=-1931;
a5=228;
Ft1=100*(a1*t^5+a2*t^4+a3*t^3+a4*t^2+a5*t);
Ft2=-500;
syms n 
a0=1/T * (int(Ft1, t, 0, 1)+int(Ft2, t, 1, 1.5));
a0 = simplify(a0);
an=2/T * (int(Ft1*cos(w*n*t),t,0,1)+ int(Ft2*cos(w*n*t),t,1,1.5));
an = simplify(an);
bn=2/T * (int(Ft1*sin(w*n*t),t,0,1)+ int(Ft2*sin(w*n*t),t,1,1.5));
bn = simplify(bn);
N = 3;
hold on
for i = 1:N
    an_subs = subs(an,n,1:i);
    bn_subs = subs(bn,n,1:i);
    F = a0+sum(an_subs.*cos(w*(1:i)*t)+bn_subs.*sin(w*(1:i)*t));
    F_handle = matlabFunction(F);
    fun = @(t,o)[o(2);(-c*a^2*o(2)-k*a^2*o(1)-m*g*a+2*F_handle(t))/(4/3*a^2*m)];
    o0 = [0 0];
    [T,O] = ode45(fun,[0 10],o0);
    plot(T,O(:,1))
end
hold off