## Matlab Assignment 9

Newbie

Posts: 23

Joined: Mon Sep 12, 2011 7:58 pm

### Matlab Assignment 9

I'm a little lost on the Newmark integration method.

1. From what I understand the first step is to use Newton's Equations of Motion just the way we have been to find qdotdot and lambda at time t=0, so that we know everything about the system at that time.

2. Then we start incrementing time. At each time step our first guess for qdotdot and lambda is the value from the previous time step and we use these values to guess updated values of q and qdot based on Newmark's equations. Based on this guess of q, we evaluate J-tilde - a value we'll use/save for every iteration in this time step.

3.From here we can use the quasi-Newton method to get a more accurate value for qdotdot and lambda using
Psi=[M*qdotdot+[Phiq]T(current_guess_qn+1)-QA(qdot,q,t);
1/(Beta*h2)*Phi(current_guess_qn+1,t)]

4. Then we update our guesses for q and qdot, check the size of the correction factor and decide whether to go through the quais-Newton again to get better values for qdotdot and lambda (back to 3 if correction is still larger than specified limit).

5. Once we're satisfied with qdotdot and lambda, we save those values along with q and qdot computed from those values and then move on to the next time step and repeat the process starting at 2.

Global Moderator

Posts: 833

Joined: Wed Sep 03, 2008 12:24 pm

### Re: Matlab Assignment 9

F11Carrie wrote:I'm a little lost on the Newmark integration method.

1. From what I understand the first step is to use Newton's Equations of Motion just the way we have been to find qdotdot and lambda at time t=0, so that we know everything about the system at that time.

2. Then we start incrementing time. At each time step our first guess for qdotdot and lambda is the value from the previous time step and we use these values to guess updated values of q and qdot based on Newmark's equations. Based on this guess of q, we evaluate J-tilde - a value we'll use/save for every iteration in this time step.

3.From here we can use the quasi-Newton method to get a more accurate value for qdotdot and lambda using
Psi=[M*qdotdot+[Phiq]T(current_guess_qn+1)-QA(qdot,q,t);
1/(Beta*h2)*Phi(current_guess_qn+1,t)]

4. Then we update our guesses for q and qdot, check the size of the correction factor and decide whether to go through the quasi-Newton again to get better values for qdotdot and lambda (back to 3 if correction is still larger than specified limit).

5. Once we're satisfied with qdotdot and lambda, we save those values along with q and qdot computed from those values and then move on to the next time step and repeat the process starting at 2.

Carrie - Actually, this is it.
One more thing: don't look at the value of \delta lambda to understand whether you need to keep iterating or not. Rather, take a look at the infinity norm of \delta qdotdot. If it's less then 1E-3, you can bail out of the quasi-Newton.
I hope this helps.
Dan

Newbie

Posts: 29

Joined: Mon Sep 12, 2011 7:58 pm

### Re: Matlab Assignment 9

My solution for k=2 is unstable. The acceleration increases to infinity. Do I need an initial guess for acceleration and lambda? I am using zeros for now.

The order of operations in my loop are:
while norm>1*10^3
1.recalculate q(n+1) and q_d(n+1) using updated q_dd(n+1)
2.recalculate phi and the jacobian using the newest values for q(n+1) q_d(n+1) and q_dd(n+1)
3.J tilde calculated only in the first iteration (or in every iteration, I've tried both)
4.Calculate the residual of q_dd and lambda
5.Calculate the correction
6.Update q_dd(n+1) and lambda(n+1) by subtracting the correction
7.Calculate the norm of the q_dd correction
end

with lambda(1)=[0 0 0 0]'
and q_dd(1)=[0 0 0 0]'
after the first iteration, my residual is
[0 19.62 -2.5 -5 19.62 0 0 0 0 0]'
which is -Qa and zeros for the scaled phi
after the first iteration, my correction is
[0 39.24 -1.6675 -10 39.24 0 0 22.12 5 -2.5]'

For the Newmark method my beta is 0.3025 and gamma is 0.6

Global Moderator

Posts: 833

Joined: Wed Sep 03, 2008 12:24 pm

### Re: Matlab Assignment 9

F11Jonathan wrote:My solution for k=2 is unstable. The acceleration increases to infinity. Do I need an initial guess for acceleration and lambda? I am using zeros for now.

The order of operations in my loop are:
while norm>1*10^3
1.recalculate q(n+1) and q_d(n+1) using updated q_dd(n+1)
2.recalculate phi and the jacobian using the newest values for q(n+1) q_d(n+1) and q_dd(n+1)
3.J tilde calculated only in the first iteration (or in every iteration, I've tried both)
4.Calculate the residual of q_dd and lambda
5.Calculate the correction
6.Update q_dd(n+1) and lambda(n+1) by subtracting the correction
7.Calculate the norm of the q_dd correction
end

with lambda(1)=[0 0 0 0]'
and q_dd(1)=[0 0 0 0]'
after the first iteration, my residual is
[0 19.62 -2.5 -5 19.62 0 0 0 0 0]'
which is -Qa and zeros for the scaled phi
after the first iteration, my correction is
[0 39.24 -1.6675 -10 39.24 0 0 22.12 5 -2.5]'

For the Newmark method my beta is 0.3025 and gamma is 0.6

Use the values from the previous time step to seed qdotdot and lambda. At time t=0 solve the EOM in conjunction with the Acc Kin Constr Eq to get qdotdot_0 and \lambda_0.
Also, in "4.Calculate the residual of q_dd and lambda", the residual is the value of Psi, not clear exactly what you mean by "residual of q_dd and lambda".
Make sure your step-size h is small enough.
How many constraints and how many generalized coordinates do you have? Do you pose the problem right?
Dan

Newbie

Posts: 29

Joined: Mon Sep 12, 2011 7:58 pm

### Re: Matlab Assignment 9

It's working now

At time t=0 solve the EOM in conjunction with the Acc Kin Constr Eq to get qdotdot_0 and \lambda_0

I was just using zeros for the initial guess

Newbie

Posts: 26

Joined: Mon Sep 12, 2011 7:58 pm

### Re: Matlab Assignment 9

Question 1
#1 Plots.jpg (38.08 KiB) Viewed 12437 times

Newbie

Posts: 22

Joined: Mon Sep 12, 2011 7:58 pm

### Re: Matlab Assignment 9

I'm stuck on how to use newton-raphson in conjugation with the backwards Euler method in problem 2. Once we run our initial guess through the equation, how do we calculate the correction factor? I thought we just divide our answer by the derivative of the equation, but this is giving me messy results. Any pointers?

Global Moderator

Posts: 833

Joined: Wed Sep 03, 2008 12:24 pm

### Re: Matlab Assignment 9

F11Christopher wrote:I'm stuck on how to use newton-raphson in conjugation with the backwards Euler method in problem 2. Once we run our initial guess through the equation, how do we calculate the correction factor? I thought we just divide our answer by the derivative of the equation, but this is giving me messy results. Any pointers?

Let's say that you are at t0 and want to find y1 at t1. Then you discretize your ODE and end up with an equation that looks like \Phi(y1)=y1 - h*f(t1,y1)=0.
Now you use Newton Raphson to find y1.
a) Take an initial guess (value of y at previous time step should be good) for y1 and set t = t+h.
b) Compute the residual corresponding to the current value of y1
c) Compute the Jacobian (sensitivity) \Psi_y1 = 1- h*partial f/partial y
d) Compute the correction \delta y = residual/Psi_y1
e) Apply the correction \delta y: y_improved = y - \delta y.
f) If correction small enough move on to the next time step (go to a) above). If not, go to b) above.

I think that the most important step is to understand the form of the nonlinear equation \Phi(y1)=y1 - h*f(t1,y1)=0 and where it comes from. After that you simply apply Newton-Raphson.
I hope this helps.
Dan

Newbie

Posts: 29

Joined: Mon Sep 12, 2011 7:58 pm

### Re: Matlab Assignment 9

Here is what I have so far.
Problem 1 Solution
Problem1.jpg (15.66 KiB) Viewed 12394 times

Problem 1 between solution and ODE 45 solution
Problem1Bonus.jpg (15.17 KiB) Viewed 12394 times

Problem 2 Forward Euler and Backward Euler Solutions
Problem2-1.jpg (14.54 KiB) Viewed 12394 times

Newbie

Posts: 29

Joined: Mon Sep 12, 2011 7:58 pm

### Re: Matlab Assignment 9

Problem 2 Difference Between Forward and Backward Euler Solutions
Problem2-2.jpg (15.03 KiB) Viewed 12394 times

Problem 3 Order Analysis
Problem3.jpg (13.9 KiB) Viewed 12394 times

Problem 4 BDF Solution
Problem4.jpg (13.79 KiB) Viewed 12394 times

Newbie

Posts: 23

Joined: Mon Sep 12, 2011 7:58 pm

### Re: Matlab Assignment 9

My solutions for the Forward and Backward Euler Methods
Problems 1 and 2.jpg (93.16 KiB) Viewed 12388 times

Problem 3 - Order analysis of Forward Euler vs. Runge-Kutta
Problem_3.jpg (76.87 KiB) Viewed 12388 times

Problem 4 - 2nd Order BDF
Problem_4.jpg (21.21 KiB) Viewed 12388 times

Newbie

Posts: 23

Joined: Mon Sep 12, 2011 7:58 pm

### Re: Matlab Assignment 9

Problem 5 - Newmark integration results
Problem_5.jpg (132.8 KiB) Viewed 12387 times

Note that torque produced by the revolute joint at the location of the joint is zero because the reaction forces pass through the point where torque is calculated. This value would be non-zero if the torque was calculated at the centroid.

Newbie

Posts: 29

Joined: Mon Sep 12, 2011 7:58 pm

### Re: Matlab Assignment 9

Here is Problem 5 Force and Torque
Problem5-F.jpg (19.52 KiB) Viewed 12382 times

Problem5-T.jpg (13.29 KiB) Viewed 12382 times

It looks like they agree with Carrie's

Newbie

Posts: 26

Joined: Mon Sep 12, 2011 7:58 pm

### Re: Matlab Assignment 9

#1 Plots.jpg (38.08 KiB) Viewed 12371 times

#2 Plots.jpg (37.39 KiB) Viewed 12371 times

Prob 3.jpg (36.74 KiB) Viewed 12371 times

Newbie

Posts: 26

Joined: Mon Sep 12, 2011 7:58 pm

### Re: Matlab Assignment 9

#4 Plot.jpg (18.83 KiB) Viewed 12371 times

Prob 5.jpg (50.68 KiB) Viewed 12371 times
Next