HW_9/Code_1.m
%% Problem no.1 (b)
% clear screen
clc
clear all
%% Newtons's Method to solve non linear equation
% initial guess
x0=[0.9;2.5;1.35;1.6];
% tol
tol = 10^-10;
% maximum iteration
nmax = 200;
% subintervals
N=100;
% anomly. function
J =@(x) gauss_J(x);
F =@(x) gauss_F(x,N);
x1 = x0;
n=0;
delF=1;
while tol x2 = x1 - inv(J(x1)) * F(x1);
x1=x2;
n=n+1;
end
for i=1:4
fprintf('x%d = %10.3f\n',i,x1(i))
end
fprintf('No. of iteration: %d\n',n)
HW_9/Code_2.m
%% Problem no. 2
% clear screen
clc
clear all
%% piece 1 clamped boundary conditions 1 to -0.67
t1 = [1 2 5 6 7 8 10 13 17];
y1 = [3 3.7 3.9 4.2 5.7 6.6 7.1 6.7 4.5];
cs1 = spline(t1,[1 y1 -0.67]);
xx = linspace(1,17,50);
plot(t1,y1,'or',xx,ppval(cs1,xx),'-');
hold on
%% piece 2 clamped boundary conditions 3 to -4.0
t2 = [17 20 23 24 25 27 27.7];
y2 = [4.5 7.0 6.1 5.6 5.8 5.2 4.1];
cs2 = spline(t2,[3 y2 -4]);
xx = linspace(17,27.7,50);
plot(t2,y2,'or',xx,ppval(cs2,xx),'-');
hold on
%% piece 3 clamped boundary conditions 0.33 to -1.5
t3 = [27.7 28 29 30];
y3 = [4.1 4.3 4.1 3.0];
cs3 = spline(t3,[0.33 y3 -1.5]);
xx = linspace(27.7,30,50);
plot(t3,y3,'or',xx,ppval(cs3,xx),'-');
axis([0,30,0,30])
grid on
%% f(x) by the c-CS
disp('c-CS')
fprintf('\nat x = 10.5, f(x) is %5.2f\n',ppval(cs1,10.5))
fprintf('\nat x = 20.5, f(x) is %5.2f\n',ppval(cs2,20.5))
fprintf('\nat x = 28.5, f(x) is %5.2f\n',ppval(cs3,28.5))
%% cubic spline natural boundary condition
pp1 = csape(t1,y1,'variational');
figure;
fnplt(pp1)
hold on
pp2 = csape(t2,y2,'variational');
fnplt(pp2)
pp3 = csape(t3,y3,'variational');
fnplt(pp3)
%% f(x) by the n-CS
disp('n-CS')
fprintf('\nat x = 10.5, f(x) is %5.2f\n',ppval(pp1,10.5))
fprintf('\nat x = 20.5, f(x) is %5.2f\n',ppval(pp2,20.5))
fprintf('\nat x = 28.5, f(x) is %5.2f\n',ppval(pp3,28.5))
%% part (b)
% lenght L for the 2nd piese of graph
ln2=0;
for i=1:numel(t2)-1
ln2 = sqrt((t2(i+1)-t2(i)).^2 + (y2(i+1)-y2(i)).^2)+ln2;
end
fprintf('\nLength L for 2nd piece : %0.3f\n',ln2)
%% part (c)
pt=linspace(17,27.7,4);
% mark the points
plot(pt(2),7.08,'.r','markersize',20)
plot(pt(3),5.59,'.r','markersize',20)
fprintf('x1 = %0.2f\n',pt(2))
fprintf('x2 = %0.2f\n',pt(3))
HW_9/Code_3.m
%% Problem 3. Halley's comet trajectory
% clear screen and close figure
clc
clear all
close all
%% solving IVP
tSpan=[0 100];
y0=[0.325514 -9.096111 -0.459460 -6.916686 0.166229 -1.305721];
options = odeset('RelTol',1e-8,'AbsTol',1e-8);
[tSol,ySol]=ode45(@cometa,tSpan,y0,options);
figure;
plot3(0,0,0,'*r','markersize',14)
hold on
plot3(ySol(:,1),ySol(:,3),ySol(:,5),'*b','markersize',2)
xlabel('X')
ylabel('Y')
zlabel('Z')
title('Trajectory Halley''s Comet')
view(-97.92,23.64)
%% t for halley's comet close to sun
rdi = sqrt(ySol(:,1).^2+ySol(:,3).^2+ySol(:,5).^2);
[rdi,Idx] = min(rdi);
disp('position of comet the time t*')
fprintf('x = %0.2f, y = %0.2f & z = %0.2f\n',ySol(Idx,1),ySol(Idx,3),ySol(Idx,5))
%% compute area between curve y and z tE(70,80)
Areay=trapz(ySol([830:1:948],3));
Areaz=trapz(ySol([830:1:948],5));
fprintf('\nArea between the curve : %0.3f\n',abs(Areay-Areaz))
%% function for calculating trajectory of comet
function dYdt=cometa(~,Y)
r=sqrt((Y(1).^2+Y(3).^2+Y(5).^2));
dxdt=Y(2);
dvxdt=-4*pi^2.*Y(1)./r.^3;
dydt=Y(4);
dvydt=-4*pi^2.*Y(3)./r.^3;
dzdt=Y(6);
dvzdt=-4*pi^2.*Y(5)./r.^3;
dYdt=[dxdt;dvxdt;dydt;dvydt;dzdt;dvzdt];
end
HW_9/gauss_F.m
%% Problem no.1 part(a)
function F=gauss_F(x,n)
% limits
a=x(1);
b=x(2);
% range
lim=linspace(a,b,n);
% eq 1
I1=0;
for i=1:numel(lim)-1
f1=@(t) sin(x(3)*(0.5*(lim(i+1)+lim(i))+0.5*(lim(i+1)-lim(i))*t)^2+x(4));
I1 = 0.5*(lim(i+1)-lim(i))*(f1(-1/sqrt(3))+f1(1/sqrt(3)))+I1;
end
I1=I1+2.886389616512186e-1;
% eq 2
I2=0;
for i=1:numel(lim)-1
f2=@(t) cos(x(3)*(0.5*(lim(i+1)+lim(i))+0.5*(lim(i+1)-lim(i))*t)^2+x(4));
I2 = 0.5*(lim(i+1)-lim(i))*(f2(-1/sqrt(3))+f2(1/sqrt(3)))+I2;
end
I2=I2-3.0709299414495893e-1;
% eq 3
I3=0;
for i=1:numel(lim)-1
f3=@(t) t*cos(x(3)*(0.5*(lim(i+1)+lim(i))+0.5*(lim(i+1)-lim(i))*t)^2+x(4));
I3 = 0.5*(lim(i+1)-lim(i))*(f3(-1/sqrt(3))+f3(1/sqrt(3)))+I3;
end
I3=I3-6.104465468269571e-1;
% eq 4
I4 = sin(x(3)*x(2)+x(4))+9.76416010290649e-1;
% Function
F=[I1;I2;I3;I4];
end
HW_9/gauss_J.m
%% Problem No. 1 (a)
function J=gauss_J(x)
% jacobian of non linear function
J = [(x(1)/2 - x(2)/2)*(x(3)*cos(x(3)*(x(1)/2 + x(2)/2 - (3^(1/2)*(x(1)/2 - x(2)/2))/3))*(3^(1/2)/6 - 1/2) - x(3)*cos(x(3)*(x(1)/2 + x(2)/2 + (3^(1/2)*(x(1)/2 - x(2)/2))/3))*(3^(1/2)/6 + 1/2)) - sin(x(3)*(x(1)/2 + x(2)/2...