mx''=-kx+Gsin(gx)e_1
e_1 1st unit vector Take m=k then the solution Is... EVeryone should know it... Anyway in case g not 1 (or in general not sqrt(k/m) we have bounded solution otherwise not (why?) we have a resonance which is bad..the floors shifts more and more and the bldg collapses.. Therefor we want bldgs with freqs out of the range of the freq of standard earthquake
int ip[3], float x[6];
ip[2]=-12;x[5]=12.34;
struct complex {
float re;
float im;};
struct complex z;
z.re=0.0;z.im=1.0;
int *k,i=0;
..
k=&i;*k=23;
printf("%d\n",i);
for(k=0;k<10;k++){
...
}
while(..){
....
} ;
do ... while(...) etc
if(...)
{ ..}
else
{...}
==, !=, ! as NOT etc
int k=(int)3.4;
Hereint k=3.4;
would give the same result but the 1st version makes the code clearer.double square(double x){ return(x*x);}
is a function and#define square1(x) (x)*(x)
is a macro computing the square of x... what is the difference? Hint: can we apply both the function and the macro to any numerical type e.g. float or integer one?fopen(),fclose(), fscanf(),fprintf(), fgetchar(), fread(),fwrite(), feof() etc
-u''+b(t)u'+c(t)u=f(t) with u(a)=al,u(b)=be
(a model linear BVP in 1d - Dirichlet boundary conditions)u(t,s)=u0(t)+s*u1(s)
substitutitng at t=b we getu(b,0)=u0(b) and u(b,1)=u0(b)+u1(b)
i.e.u1(b)=u(b,1)-u(b,0).
Now we get s0 solving a liner scalar equation:be=u(b,s0)=u0(b)+s0*u1(b)
gettings0=(be-u0(b))/u1(b)=(be-u(b,0))/(u(b,1)-u(b,0)).
Naturally, u(b,0) and u(b,1) are computed by some numerical approximated scheme. If computed u1(b) is nonzro (in practise above some level; in fl we almost always get nonzero due to rounding errors, an ODE solver add some error..), then we may be almost sure that the problem has a unique solution (almost as our ODE solver may always fail miserably).-u''+c(x)u=f
with Dirichelt left bnd cnd:u(a)=al
and the right Neumann cond:u'(b)=be
(-u_{N-1}+2u_n-u_{N+1})/h^2+c(b)u_N=f(b)
(u(b)-u(b-h))/h-be=-u''(b)h/2 + O(h^2)
but then using our equation:-u''(b)=f(b)-c(b)u(b)
so using this we get that(u(b)-u(b-h))/h -be= (f(b)-c(b)u(b))h/2 + O(h^2)
so replacing u(b) by u_n and u(b-h) by u_n we get a second order scheme:(u_N-u_{N-1}))/h +0.5c(b)u_N*h=be + f(b)*h/2
(please check if everything is OK...)u(b)=be
We then consider a new space Vl with functions being zero at the left end (Dirichlet bnd condition - in general in FEM- Dirichlet conditions go into the definition of the space) and the weak formulation (again integrating by parts):u \in Vl \int_a^b u'v' +c uvdx = \int_a^b f v dx +u'(b)v(b) for all v in Vl
New linear FEM space (functions that are zero at the left end only):Vlh=\{u in C: u(a)=0,u linear on [x_k,x_{k+1}]\}
and variational problem:u \in Vlh \int_a^b uh'v' +c uhvdx = \int_a^b f v dx +be*v(b) for all v in Vlh
Nodal basis (gk) k=1,...,N (one more extra nodal function relatedto the end b!) and we get the systemAl * U =Fl
Al has an extra row and column (but is still symmetric and positive definite and tridiagonal - problem) Fl has also an extra entry i.e Fl=[F;Fl(N)] the last one namely: Fl(N)= \int_a^b f g_N dx + be*g_N(b) which can be approximated (integral by trapezoidal rule on [x_{N-1},x_N]) as 0.5f(b)(x_N-x_{N-1})+be in case of equaidistant mesh: 0.5*f(b)h +be - we get that this equation (the whole one not just the rhs..) is equivalent to the FDM approximation equation of Neumann cond by the backward difference + correction by utilizing the ODE equation at b...s=0; for(k=0;k<16;k++) s+=x[k];
one can write
s=x[0]+...+x[15];
for(k=0;k<16;k++) x[k]=x[k]/a;
one can write
ia=1./a;(k=0;k<16;k++) x[k]=ia*x[k];
G*s=r
with G=[X'*X,X'*Z;Z'*X,k*A^(-1)] and r=[X'*y;Z'*y]y''=-y with y(0)=1;y'(0)=0 on [0,2*pi](solution: a*cos(x)+bsin(x))
y''=-sin(y)
with the same initial data. x'=998*x+1998*y;
y'=-999*x-1999*y ;
x(0)=y(0)=1 on [0,10], [0,100] etc
[998 1998;-999, -1999]
x'=-0.04x+10^4*y*z;
y'=0.04*x-10^4*y*z-3*10^7y*y;
z'=3*10^7*y*y
y'=200y with y(0)=1 on [0,10]
What happens? Repeat fory'=y*y with y(0)=1 on [0,10]
i.e. with IVP with a real blowup.. Think if we can distinguish between them? I.e. If we get blowup we wouldike to know if it is a real one or a numerical one.y''=-y with y(0)=1;y'(0)=0 on [0,T]
for T=2k\Pi k=1,2,10,100 etc check if |y(T)-1|==0.y''=-sin(y) with y(0)=1 on [0,100]
The solution is also periodic but with unknown period.. How to find it numerically? (ok it is possible to compute it analytically...) y'=f(y,t),
y(0)=s
y'=y*(1-y)*sin(y*y+x);
y(0)=s.
m_i*x_i''=k_{i+1}(x_{i+1}-x_i) -k_i(x_{i}-x_{i-1}), forces on the i-th floor, i>0
on the zero floor (ground floor) acts also analogous restorative force proportional to the shift relative to the ground: i.e.-k_0x_0
and the restorative force between zero floor and 1st floor: i.e.m_0x_0''=k_1(x_1-x_0) -k_0x_0.
Naturally, on the last floor n-th acts only one restorative force...m_n*x_n''=-k_n(x_n-x_{n-1}).
f(t)=G*sin(gamma t)
acts on the ground floor (level 0) e.g .gamma=3.
MX'=SX+F
for M diagonal matrix n+1 x n+1 (masses of floors on the diagonal), matrix S a 3-diagonal symmetric matrix such that on the main diagonal we have -k_{i+1}-k_i (k_{n+1}=0) sub- super- diagonal = (k_1,..,k_n)for(){ }, while(){ }, do{ } while()
,
option - p should be given interactively from the keyboard by an user
fopen(), fclose(), fwrite(), fread(),
fprintf(), fscanf(), fgetchar() feof()etc
Write a code that will write a vector (double
array)to a file - binary or ascii (text) file
then another code that will read it from the file.
Test for: [1 2 3 4 -9.45].
float,double
etc.
Write a a C function generating a vector sin(x)
for a given vector x (e.g. equidistant points on [a,b]) -
parameters : a pointer to x (or ends of interval - easier version))
and pointer to returned vector with sin(x) - C function should allocate memory
and write to that memory sin(x), optionally export this vector to
binary file (fwrite()
)
and then read it back to memory to another memory address.
#define #if
etc
Modify the code from the previous function in such a way that depending on the precompiler options it will work either for doubles or floats (single
precision).
#include
"myblas.h"
) with the proper headers of this functions (in order to use them in your C code)
float
array in C. Then show the matrix on the computer display, i.e. save a matrix in text format in octave and then display by your C programm - check if it is OK.
(2)Write another C code which writes a matrix from C to text file in a default text format of octave. Read it using octave - check if it works.
(3)Try another wrtitng to other octave formats (ascii files) - or to binary octave format - (difficult - one has to find details of octave binary format)
dgtsv()
solve the tridiagonal system
system Ax=f with the tridiagonal matrix A with 7 on the main diagonal, -2 on superdiagonal,
and -1 on subdiagonal for f=[7;-1;0;0;..;0] for N=15.Check if you get the right solution!
dnrm2()
and ddot()
dnrm2()
function - and dynamically allocating memory (using
malloc(); free()
functions
dgesv()
function - solving a simple linear system and computing the inverse of A
dgels()
function solving a regular LLSP with i.e.
max rank matrix A MxN - underdetermined system i.e. M>N; there is also a function linreg()
solving linear regression using dgels()
-u''+u=0; u(0)=u(T)=1
Apply the shooting method, namely, solve instead an Initial Value Problem (IVP)-u''+u=0; u(0)=1; u'(0)=s
getting u(T,s) since it depends on s, (practically numerically using lsode()). Find s0 s.t. u(t;s0)=1. It is enough to shoot twice i.e. take s2=1 and s1=0 and compute using lsode() u(t;0) and u(t;1). Then our s0=(1-u(T,0))/(u(T,1)-u(T,0)) (why?). Then compute u(T;s0) using lsode() and |u(T;s0)-1|. Why is the error growing with T?u''=(3/2)u^2 with u(0)=4 and u(1)=1.
Define using lsode() the function F(s)=u(1;s)-1 with u(t;s) being the solution of IVP -u''=(3/2)u^2 with u(0)=4 and u'(0)=s. Solve the system F(s*)=0 using fzero() of fsolve() taking different starting s0 initial iterates for a nonlinear solver Plot the graph of F on [-100,0].-u''+c*u =f
u(a)=alpha, u(b)=beta
uh(k) approximating u(x(k))
with x(k)=a+k*h h=(b-a)/N - uh(0)=alpha uh(N)=beta - solve the linear systemAU=F with U=(uh(0),..,uh(N)) F=[alpha,f(x(1)),..,f(x(N-1)),beta]
andA=A1+c*Id
- A1 tridiagonal with A1(0,0)=A1(N,N)=1, A1(i,i)=2/h^2 and A1(i-1,i)=A1(i+1,i)=-1/h^2 i=1,..,N-1.-u''=f
u(0)=al;u'(1)=be;
(u_{N-1}-u_{N+1})/(2h)=be
the remaining FDM equations are:
u_h(x_0)=al
(-u_{k-1}+2u_k-u_{k+1})/(h^2)=f(x_k) k=1,...,N+1
-y''=1 y(0)=s0;y'(1)=0;
H- given function on [0.5,1], (simplified project i.e. 1d instead 2d) Solve the problem by FDM method on the equidistant mesh x(k)=k*h k=0,..,N with h=1/N.u_t=u_{xx} in (0,T] x (0,L)
- space 1D -right bnd homogeneous Neumann condition:u'(t,L)=0
left bnd condition : Dirichlet -u(t,0)=s(t)
the initial condition:u(0,x)=u0(x) given
. Solve by FDM using lsode() after space FDM dicretization or solve ODE by using implicit Euler scheme.u(t,x(k))'- (u(t,x(k-1))+2(u(t,x(k))- (u(t,x(k+1)))/(h^2)=f(t,x(k)) k=1,..,N-1
Discrete FDM left Dirichlet bnd condition:u(t,x(0))=s(t);
Discrete FDM right Neumann zero boundary codition:[u(t,x(N))-u(t,x(N-1))]/h=0;
Discrete initial condition:u(0,x(k))=u0(x(k)) k=1,..,N-1;
Mesh:x(k)=k*h h=L/N;
f(s)=\int_0^T \int_c^b|y(t,x;s)-H(t,x)|^2 dx
is minimal for s0 where y solves 1D parabolic PDE:y_t-y_{xx}=0 on (0,b)
with the given initial condition:y(0,x)=y0(x) given,
and boundary values:y_x(t,b)=0 and y(t,0)=s0(t) (control),
Control is s0 - i.e .we can control only the left Dirichlet boundary condition in t.H- given function on [c,b] (0 < c < b),
dnrm2()
and ddot()
dnrm2()
function - and dynamically allocating memory (using
malloc(); free()
functions
dgesv()
function - solving a simple linear system and computing the inverse of A
dgels()
function solving a regular LLSP with i.e.
max rank matrix A MxN - underdetermined system i.e. M>N; there is also a function linreg()
solving linear regression using dgels()